C
C *** LAST REVISED ON 12-AUG-1987 11:58:37.75
C *** SOURCE FILE: [DL.GRAPHICS.LONGLIB]HIDELIB.FOR
C
C *************************************************************************
C
C	ROUTINES FOR 3D PLOTTING COMPATIBLE WITH HIDDEN LINE REMOVAL
C	ROUTINES FOR LONGLIB GRAPHICS LIBRARY.
C	THIS FILE CONTAINS SOME SUBROUTINES ORIGINALLY FROM COSMIC THAT HAVE
C	BEEN MODIFIED TO BE COMPATIBLE WITH THE LONGLIB GRAPHICS LIBRARY.
C
C	ROUTINES IN THIS FILE ARE FORTRAN 77 COMPATIBLE WITH THE
C	FOLLOWING EXCEPTIONS:
C		1. TABS (^I) ARE USED TO INDENT LINES
C		2. TRAILING COMMENTS SEPARATED BY ! ARE USED
C		3. BYTE VARIABLE TYPES ARE USED IN NUM3DH AS ARE
C		   ENCODE
C		4. INTEGER*2 IS USED IN SYM3DH TO SAVE STORAGE SPACE
C
C ************************************************************************
C
	SUBROUTINE INIT3DH(X,Y,Z,YA,RO,PI,ZOM,DIS,MNE,MNH)
C
C	CREATED: DGL 1-AUG-1984
C
C	INITALIZE 3D HIDDEN ROUTINES
C
C	X,Y,Z	(R)	RELATIVE ORIGIN OFFSET
C	YA,RO,PI(R)	YAW, ROLL, PITCH EULERIAN ANGLES IN DEGREES
C	ZOM	(R)	RELATIVE SCALE FACTOR
C	DIS	(R)	VIEWING DISTANCE  9999. = INFINITY (NO PERSPECTIVE)
C	MNE,MNH (I)	MAXIMUM NUMBER OF EDGES AND HOLES ON ONE POLYGON
C
	COMMON /CHIDE/AL(3),V0(3),VSF,SCX,YAW,ROL,PIT,LZ,VP,JJJ
	V0(1)=X
	V0(2)=Y
	V0(3)=Z
	YAW=YA		! PSI
	ROL=RO		! PHI
	PIT=PI		! THETA
	VSF=ABS(ZOM)
	VP=ABS(DIS)
	SCX=SIGN(1.,ZOM)! HIDDEN LINE PLOT UNITS PER INCH
C			! IF ZOM < 0, ALL LINES SHOWN
	LZ=IABS(MNE)	! MAX NUMBER OF EDGES ON ONE POLYGON
	IF (LZ.LT.4) LZ=4
	JJJ=LZ+2+2*MNH	! MNH IS MAX NUMBER OF HOLES IN ONE POLYGON
	RETURN
	END
C
	SUBROUTINE PLT3DH(X,Y,Z,IP)
C
C	CREATED: DGL 1-AUG-1984
C
C	3D PLOT ROUTINE COMPATABLE WITH HIDDEN LINE REMOVAL ALGORITHM
C
C	X,Y,Z	(R)	POINT
C	IP	(R)	CONTROL FLAG
C			= 0	CHANGE COLOR TO X, RELATIVE PLOTTING ANGLE TO Y
C				CLEAR SCREEN WHEN X < 0
C				CHANGE RELATIVE SCALE FACTOR BY Z
C				(RESET SCALE FACTOR TO 1. WHEN Z < 0.)
C				(S.F. NOT CHANGED WHEN Z=0)
C			= 2	PLOT TO (X,Y,Z) PEN DOWN
C			=-2	PLOT TO (X,Y,Z) PEN DOWN, NEW ORIGIN
C			= 3	MOVE TO (X,Y,Z) PEN UP
C			=-3	MOVE TO (X,Y,Z) PEN UP, NEW ORIGIN
C			= 9	ERASE TO (X,Y,Z) PEN UP
C			=-9	ERASE TO (X,Y,Z) PEN UP, NEW ORIGIN
C
	COMMON /CHIDE/AL(3),V0(3),VSF,SCX,YAW,ROL,PIT,LZ,VP,JJJ
	EQUIVALENCE (AL(1),XJ),(AL(2),YJ),(AL(3),ZJ)
	DATA PI/0.0174533/
C
C	STORE EULERIAN ANGLES.
C
	XX=YAW*PI
	YY=ROL*PI
	ZZ=PIT*PI
	COSY=COS(YY)
	SINY=SIN(YY)
	COSZ=COS(ZZ)
	SINZ=SIN(ZZ)
	COSX=COS(XX)
	SINX=SIN(XX)
C
C	TRANSFORM (X,Y,Z)
C
	XJ=VSF*X+V0(1)
	YJ=VSF*Y+V0(2)
	ZJ=VSF*Z+V0(3)
C
	X1=ZJ*(COSY*SINX)+XJ*(COSY*COSX)-YJ*SINY
	TW=YJ*COSY*COSZ
	TZ=XJ*(SINZ*SINX+SINY*COSZ*COSX)
	TY=ZJ*(-SINZ*COSX+SINY*COSZ*SINX)
	Y1=TZ+TW+TY
	PT=YJ*COSY*SINZ
	PK=ZJ*(COSZ*COSX+SINY*SINZ*SINX)
	PS=XJ*(-COSZ*SINX+SINY*SINZ*COSX)
	ZJ=PK+PS+PT
	VL=ABS(VP)
C
CC	IF (VP.EQ.9999.) GOTO 5
C
C     CALCULATES PERSPECTIVE BASED ON VALUE VP(DV) 
C
	HH=VL/(VL-ZJ)
	XJ=X1*HH
	YJ=Y1*HH
	ZJ=ZJ*HH
5	CONTINUE
C
	IF (IABS(IP).GT.9) GOTO 8
	GOTO (10,8,20,20,8,8,8,8,8,20) IABS(IP)+1
8	RETURN
10	CONTINUE	! COLOR CHANE
	IF (Z.NE.0.0) VSF=VSF*Z
	IF (Z.LT.0.) VSF=1.0
	CALL PLOT(X,Y,IP)
	GOTO 8
20	CONTINUE	! PLOT
	IF (IP.GT.0) GOTO 25
	V0(1)=X1
	V0(2)=Y1
	V0(3)=ZJ
25	CALL PLOT(XJ,YJ,IABS(IP))
	GOTO 8
	END
C
C
	SUBROUTINE WHERE3H(X,Y)
C
C	CREATED: DGL 1-AUG-1984
C
C	RETURNS THE SCREEN COORDINATES OF THE LAST PLOTTED POINT
C	FROM PLT3DH
	COMMON /CHIDE/AL(3),V0(3),VSF,SCX,YAW,ROL,PIT,LZ,VP,JJJ
	X=V0(1)
	Y=V0(2)
	RETURN
	END
C
C
	SUBROUTINE XFRM3D(X,Y,Z,X2,Y2,Z2)
C
C	CREATED: DGL 1-AUG-1984
C
C	3D TRANSFORM (X,Y,Z) TO (X2,Y2,Z2)
C
C	X,Y,Z	(R)	INPUT POINT
C	X2,Y2,Z2(R)	OUTPUT TRANSFORMED POINT
C
	COMMON /CHIDE/AL(3),V0(3),VSF,SCX,YAW,ROL,PIT,LZ,VP,JJJ
	EQUIVALENCE (AL(1),XJ),(AL(2),YJ),(AL(3),ZJ)
	DATA PI/0.0174533/
C
C	STORE EULERIAN ANGLES.
C
	XX=YAW*PI
	YY=ROL*PI
	ZZ=PIT*PI
	COSY=COS(YY)
	SINY=SIN(YY)
	COSZ=COS(ZZ)
	SINZ=SIN(ZZ)
	COSX=COS(XX)
	SINX=SIN(XX)
C
C	TRANSFORM (X,Y,Z)
C
	XJ=VSF*X+V0(1)
	YJ=VSF*Y+V0(2)
	ZJ=VSF*Z+V0(3)
	X1=ZJ*(COSY*SINX)+XJ*(COSY*COSX)-YJ*SINY
	TW=YJ*COSY*COSZ
	TZ=XJ*(SINZ*SINX+SINY*COSZ*COSX)
	TY=ZJ*(-SINZ*COSX+SINY*COSZ*SINX)
	Y1=TZ+TW+TY
	PT=YJ*COSY*SINZ
	PK=ZJ*(COSZ*COSX+SINY*SINZ*SINX)
	PS=XJ*(-COSZ*SINX+SINY*SINZ*COSX)
	ZJ=PK+PS+PT
	VL=ABS(VP)
CC	IF (VP.EQ.9999.) GOTO 5
C
C     CALCULATES PERSPECTIVE BASED ON VALUE VP(DV) 
C
	HH=VL/(VL-ZJ)
	XJ=X1*HH
	YJ=Y1*HH
	ZJ=ZJ*HH
5	CONTINUE
	X2=XJ
	Y2=YJ
	Z2=ZJ
	RETURN
	END
C
C ************************************************************************
C
C	AUXILARY ROUTINES FOR THE INIT3DH PORTION OF THE LONGLIB
C	GRAPHICS LIBRARY
C
	SUBROUTINE AXIS3DH(X0,Y0,Z0,PA,PB,PG,T,N0,S0,B0,C0,D0,E0,F0,IOP)
C
C	CREATED: DGL FEB. 1986
C
C	X0	X COORDINATE OF START OF AXIS [DEFAULT 0.0]
C	Y0	Y COORDINATE OF START OF AXIS [DEFAULT 0.0]
C	Z0	Z COORDINATE OF START OF AXIS [DEFAULT 0.0]
C	T	CHARACTER STRING TO DESCRIBE AXIS [DEFAULT NONE]
C	PA	ANGLE FROM X-Y PLANE (DEG) [DEFAULT 0.0]
C	PB	ANGLE FROM X-Z PLANE (DEG) [DEFAULT 0.0]
C	PG	ANGLE OF ROTATION AROUND PA,PB RAY (DEG) [DEFAULT 0.0]
C	N0	NUMBER OF CHARACTERS IN STRING [DEFAULT +0]
C		- ON CLOCKWISE SIDE OF AXIS (NORMAL FOR X)
C		+ ON COUNTER CLOCKWISE SIDE OF AXIS (NORMAL FOR Y)
C		HUNDREDS DIGIT	= 1 NO LABELING OF AXIS--TICKS AND LINE ONLY
C		THOUSANDS DIGIT = 1 TICKS PLACED ON OPPOSITE SIDE OF AXIS LABEL
C	S0	LENGTH OF AXIS IN INCHES
C		< 0	AXIS RUNS BACKWARDS
C		= 0	NO ACTION (NOP)
C		> 0	NORMAL
C	B0	MINIMUM VALUE ON TICK AXIS
C	C0	MAXIMUM VALUE ON TICK AXIS
C	D0	INT(D0) = NUMBER OF MAJOR AXIS TICKS
C		INT((INT(D0)-D0)*100) = NUMBER OF MINOR AXIS TICKS BETWEEN MAJOR TICKS
C		INT(MOD((INT(D0)-D0)*10000,100) = NUMBER OF SUB MINOR AXIS TICKS 
C			BETWEEN MINOR TICKS
C	E0	CHARACTER SIZE OF TITLE AND NUMBER (IF E0=0 DEFAULTS TO .15)
C			< 0 THEN DO NOT AUTO SCALE BY (x10 TO POWER)
C	F0	NUMBER SPECIFICATION (FORMAT FX.Y) (2 DIGITS FOR Y, MAX(X)=18)
C			[DEFAULT 6.02]
C	IOP	FLAG TO INDICATE SKETCH OR PLT3DH USAGE
C		IOP=0: PLT3DH
C		IOP=1: SKETCH
C
	DIMENSION ROT(4,4),TEMP(4,4),V(4)
	DIMENSION XX(2),YY(2),ZZ(2)
	INTEGER T(1)
	DATA SPACE/0.08/		! MIN SPACING BETWEEN ITEMS
	XROT3D(X,Y)=XVMUL3D(1,X,Y,0.0,V,ROT)+X01
	YROT3D(X,Y)=XVMUL3D(2,X,Y,0.0,V,ROT)+Y01
	ZROT3D(X,Y)=XVMUL3D(3,X,Y,0.0,V,ROT)+Z01
C
	TL=0.1				! TICK LENGTH
	X01=X0
	Y01=Y0
	Z01=Z0
	PAA=PA
	PBA=PB
	PGA=PG				! ROTATION ANGLES
	E1=E0				! CHARACTER SIZE
	IOPT=IOP
	CS=ABS(E1)
	IF (CS.EQ.0.) CS=.15
	IF (S0.EQ.0.0) GOTO 1000	! ZERO LENGTH AXIS
	X1=D0*1.000002
	NJT=ABS(X1)			! NUMBER OF MAJOR TICKS
	NNT=100.*(ABS(X1)-NJT)		! NUMBER OF MINOR TICKS
	NST=100.*((ABS(X1)-NJT)*100.-NNT)	! NUMBER OF SUB-MINOR TICKS
	IF (NJT.LT.2) NJT=2
	XJ=S0/(NJT-1)			! INCREMENT BETWEEN MAJOR TICKS
	XN=XJ
	IF (NNT.NE.0) XN=XN/(NNT+1)	! INCREMENT BETWEEN MINOR TICKS
	XS=XN/(NST+1)			! INCREMENT BETWEEN SUB-MINOR TICKS
	N1=N0
	HOR=0.0
	IF (MOD(IABS(N1),10000).GE.1000) HOR=90. ! DEFAULT ROTATION ANGLE
	NC=MOD(IABS(N1),100)		! NUMBER OF CHARACTERS IN TITLE
	FA=F0				! FORMAT NUMBER
	ND=MOD(ABS(FA),1000.)		! NUMBER OF DIGITS
	NG=100.*(MOD(ABS(FA),1000.)-ND)	! NUMBER OF DIGITS RIGHT OF D.P.
	IF (NG.GT.17) NG=NG/10		! CORRECT SIMPLE ERRORS
	IF (ND.EQ.0) ND=NG
	IF (ABS(FA).GT.1000) NG=-1	! FORMATTED INTEGER
	TL1=TL
	IF (IABS(N1)/10000.GT.0) TL1=-TL! REVERSE SIDE OF TICKS
	IF (IABS(N1)/10000.GT.0) TL=0.	! REVERSE SIDE OF TICKS
	IF (IABS(N1).GE.1000) GOTO 10
		DNX=-ND*CS/2.		! NUMBER LABEL DISTANCE FROM AXIS
		DNY=-TL-SPACE-CS	! NUMBER LABEL DISTANCE FROM AXIS
		DTY=DNY-CS-SPACE	! TITLE DISTANCE FROM AXIS
		GOTO 20
10	CONTINUE			! HORIZONTAL NUMBERS ON VERTICAL AXIS
		DNX=-CS/2.		! NUMBER LABEL DISTANCE FROM AXIS
		DNY=-TL-SPACE		! NUMBER LABEL DISTANCE FROM AXIS
		DTY=DNY-ND*CS-SPACE	! TITLE DISTANCE FROM AXIS
20	DTX=(ABS(S0)-NC*CS)/2.		! TITLE DISTANCE FROM AXIS
	IF (S0.LT.0) DTX=-DTX-CS*NC
	IF (N1.LT.0) GOTO 30		! CLOCKWISE TITLES
		DNY=-DNY-CS		! COUNTER-CLOCKWISE TITLES
		DTY=-DTY-CS
		TL1=-TL1		! CHANGE SIDES OF TICKS
		IF (IABS(N1).LT.1000) GOTO 30
			DNY=DNY+CS*ND
			DTY=DNY+SPACE
30	CONTINUE
C
C	SET UP 3-D CONVERSION MATRIXES
C
	CALL ROTEM(1,PGA,ROT)
	CALL ROTEM(3,PBA,TEMP)
	CALL MATMUL4(ROT,ROT,TEMP)
	CALL ROTEM(2,PAA,TEMP)
	CALL MATMUL4(ROT,ROT,TEMP)
C
	X1=0.0				! FIRST MAJOR TICK
	Y1=0.0
	Y2=-TL1
	IF (IOPT.EQ.0) THEN
		CALL PLT3DH(XROT3D(X1,Y1),YROT3D(X1,Y1),ZROT3D(X1,Y1),3)
		CALL PLT3DH(XROT3D(X1,Y2),YROT3D(X1,Y2),ZROT3D(X1,Y1),2)
	ELSE
		XX(1)=XROT3D(X1,Y1)
		XX(2)=XROT3D(X1,Y2)
		YY(1)=YROT3D(X1,Y1)
		YY(2)=YROT3D(X1,Y2)
		ZZ(1)=ZROT3D(X1,Y1)
		ZZ(2)=ZROT3D(X1,Y1)
		IER=0
		NN=2
		CALL SKETCH(XX,YY,ZZ,NN,IER)
	ENDIF
	DO 100 I=1,NJT-1		! MAJOR TICKS
		IF (IOPT.EQ.0) THEN
		  CALL PLT3DH(XROT3D(X1,Y1),YROT3D(X1,Y1),ZROT3D(X1,Y1),3)
		ELSE
		  XX(1)=XROT3D(X1,Y1)
		  YY(1)=YROT3D(X1,Y1)
		  ZZ(1)=ZROT3D(X1,Y1)
		ENDIF
		X1=XJ*I
		Y2=-TL1
		IF (IOPT.EQ.0) THEN
		  CALL PLT3DH(XROT3D(X1,Y1),YROT3D(X1,Y1),ZROT3D(X1,Y1),2)
		  CALL PLT3DH(XROT3D(X1,Y2),YROT3D(X1,Y2),ZROT3D(X1,Y2),2)
		ELSE
		  XX(2)=XROT3D(X1,Y1)
		  YY(2)=YROT3D(X1,Y1)
		  ZZ(2)=ZROT3D(X1,Y1)
		  IER=0
		  NN=2
		  CALL SKETCH(XX,YY,ZZ,NN,IER)
		  XX(1)=XROT3D(X1,Y2)
		  YY(1)=YROT3D(X1,Y2)
		  ZZ(1)=ZROT3D(X1,Y2)
		  IER=0
		  NN=2
		  CALL SKETCH(XX,YY,ZZ,NN,IER)
		ENDIF
		DO 110 J=1,NNT+1	! MINOR TICKS
			Y2=-TL1*.7
			X2=X1+XN*J-XJ
			IF (IOPT.EQ.0) THEN
		CALL PLT3DH(XROT3D(X2,Y1),YROT3D(X2,Y1),ZROT3D(X2,Y1),3)
		CALL PLT3DH(XROT3D(X2,Y2),YROT3D(X2,Y2),ZROT3D(X2,Y2),2)
			ELSE
		XX(1)=XROT3D(X2,Y1)
		XX(2)=XROT3D(X2,Y2)
		YY(1)=YROT3D(X2,Y1)
		YY(2)=YROT3D(X2,Y2)
		ZZ(1)=ZROT3D(X2,Y1)
		ZZ(2)=ZROT3D(X2,Y2)
		IER=0
		NN=2
		CALL SKETCH(XX,YY,ZZ,NN,IER)
			ENDIF
			Y2=-TL1*.4
			DO 120 K=1,NST	! SUB MINOR TICKS
				X3=X2+XS*K-XN
				IF (IOPT.EQ.0) THEN
		CALL PLT3DH(XROT3D(X3,Y1),YROT3D(X3,Y1),ZROT3D(X3,Y1),3)
		CALL PLT3DH(XROT3D(X3,Y2),YROT3D(X3,Y2),ZROT3D(X3,Y2),2)
				ELSE
		XX(1)=XROT3D(X3,Y1)
		XX(2)=XROT3D(X3,Y2)
		YY(1)=YROT3D(X3,Y1)
		YY(2)=YROT3D(X3,Y2)
		ZZ(1)=ZROT3D(X3,Y1)
		ZZ(2)=ZROT3D(X3,Y2)
		IER=0
		NN=2
		CALL SKETCH(XX,YY,ZZ,NN,IER)
				ENDIF
120			CONTINUE
110		CONTINUE
100	CONTINUE
	IF (MOD(IABS(N1),1000)/100.NE.0) GOTO 1000 ! NO LABELING
	XS=0.0					! EXPONENT
	IF (E1.LT.0.) GOTO 140			! NO AUTO SCALING
	X1=0.
	Y1=0.
	I=ND-NG-2
	IF (B0.LT.0..OR.C0.LT.0.) I=I-1
	IF (B0.NE.0.) X1=ALOG10(ABS(B0))
	IF (X1.LT.0..AND.X1.NE.AINT(X1)) X1=X1-1.
	X1=AINT(X1)
	IF (C0.NE.0.) Y1=ALOG10(ABS(C0))
	IF (Y1.LT.0..AND.Y1.NE.AINT(Y1)) Y1=Y1-1.
	Y1=AINT(Y1)
	X2=MIN(X1,Y1)
	X3=MAX(X1,Y1)
	IF (X2.LT.0.0.AND.NG.LE.-X2) XS=-NG+1-X2
	IF (I.LE.X3+XS) XS=XS+(I-X3-XS)-1
140	Y1=DNY
	Y2=(C0-B0)/(NJT-1)
	E1=XS					! EXPONENT VALUE
	DO 150 I=1,NJT				! LABEL MAJOR TICKS
		X1=(I-1)*XJ+DNX
		C1=(Y2*(I-1)+B0)*10.**E1
		CALL NUM3DH(XROT3D(X1,Y1),YROT3D(X1,Y1),ZROT3D(X1,Y1),
     $			PAA,PBA,PGA,CS,C1,FA,IOPT)
150	CONTINUE
	IF (E1.NE.0.0) DTX=DTX-CS*3.		! ADD EXPONENT SPACE
	X1=DTX
	Y1=DTY
	IF (NC.NE.0) CALL SYM3DH(XROT3D(X1,Y1),YROT3D(X1,Y1),ZROT3D(X1,Y1),
     $		PAA,PBA,PGA,CS,T,NC,IOPT)
	IF (E1.EQ.0.0) GOTO 1000		! NO EXPONENT
	X1=X1+NC*CS
	CALL SYM3DH(XROT3D(X1,Y1),YROT3D(X1,Y1),ZROT3D(X1,Y1),
     $		PAA,PBA,PGA,CS,'(X10)',5,IOPT)
	X1=X1+4.5*CS
	Y1=Y1+CS/2.
	CS=CS/2.
	CALL NUM3DH(XROT3D(X1,Y1),YROT3D(X1,Y1),ZROT3D(X1,Y1),
     $		PAA,PBA,PGA,CS,E1,-1,IOPT)
1000	RETURN
	END
C
C
	SUBROUTINE SYM3DH(XA,YA,ZA,PA,PB,PG,HI,ASC,NC,IOP)
C
C	CREATED: DGL 1-OCT-1984
C	MODIFIED: DGL FEB 1986 -- ADDED IOP WITH SKETCH CALL
C
C	PLOTS SYMBOLS IN 3-D USING PLT3DH
C
C	XA,YA,ZA COORDINATES OF STRING LOWER-LEFT CORNER
C		  (999.,999.,999.) CONTINUE FROM LAST POINT
C	PA,PB	SPECIFIES THE PROJECTED ANGLES (IN DEGS) RELATIVE TO THE
C		 X-Y PLANE AND THE X-Z PLANE OF THE LINE ALONG THE
C		 BASE OF THE PLOTTED NUMBER SEQUENCE
C	PG	SPECIFIES THE ROTATION OF THE PLOTTED SYMBOLS AROUND
C		THE LINE SPECIFIED BY PA,PB
C	HI	HEIGHT OF THE PLOTTED CHARACTERS
C	ASC	BYTE ARRAY CONTAINING CHARACTERS TO BE PLOTTED
C	NC	NUMBER OF CHARACTERS TO BE PLOTTED
C	IOP	FLAG TO INDICATE SKETCH OR PLT3DH USAGE
C		IOP=0: PLT3DH
C		IOP=1: SKETCH
C
      REAL ASC(1)
      INTEGER*2 SYM(128),STB(416)
      INTEGER*2 SYN  1(8)
      INTEGER*2 SYN  9(8)
      INTEGER*2 SYN 17(8)
      INTEGER*2 SYN 25(8)
      INTEGER*2 SYN 33(8)
      INTEGER*2 SYN 41(8)
      INTEGER*2 SYN 49(8)
      INTEGER*2 SYN 57(8)
      INTEGER*2 SYN 65(8)
      INTEGER*2 SYN 73(8)
      INTEGER*2 SYN 81(8)
      INTEGER*2 SYN 89(8)
      INTEGER*2 SYN 97(8)
      INTEGER*2 SYN105(8)
      INTEGER*2 SYN113(8)
      INTEGER*2 SYN121(8)
      INTEGER*2 STC  1(8)
      INTEGER*2 STC  9(8)
      INTEGER*2 STC 17(8)
      INTEGER*2 STC 25(8)
      INTEGER*2 STC 33(8)
      INTEGER*2 STC 41(8)
      INTEGER*2 STC 49(8)
      INTEGER*2 STC 57(8)
      INTEGER*2 STC 65(8)
      INTEGER*2 STC 73(8)
      INTEGER*2 STC 81(8)
      INTEGER*2 STC 89(8)
      INTEGER*2 STC 97(8)
      INTEGER*2 STC105(8)
      INTEGER*2 STC113(8)
      INTEGER*2 STC121(8)
      INTEGER*2 STC129(8)
      INTEGER*2 STC137(8)
      INTEGER*2 STC145(8)
      INTEGER*2 STC153(8)
      INTEGER*2 STC161(8)
      INTEGER*2 STC169(8)
      INTEGER*2 STC177(8)
      INTEGER*2 STC185(8)
      INTEGER*2 STC193(8)
      INTEGER*2 STC201(8)
      INTEGER*2 STC209(8)
      INTEGER*2 STC217(8)
      INTEGER*2 STC225(8)
      INTEGER*2 STC233(8)
      INTEGER*2 STC241(8)
      INTEGER*2 STC249(8)
      INTEGER*2 STC257(8)
      INTEGER*2 STC265(8)
      INTEGER*2 STC273(8)
      INTEGER*2 STC281(8)
      INTEGER*2 STC289(8)
      INTEGER*2 STC297(8)
      INTEGER*2 STC305(8)
      INTEGER*2 STC313(8)
      INTEGER*2 STC321(8)
      INTEGER*2 STC329(8)
      INTEGER*2 STC337(8)
      INTEGER*2 STC345(8)
      INTEGER*2 STC353(8)
      INTEGER*2 STC361(8)
      INTEGER*2 STC369(8)
      INTEGER*2 STC377(8)
      INTEGER*2 STC385(8)
      INTEGER*2 STC393(8)
      INTEGER*2 STC401(8)
      INTEGER*2 STC409(8)
      EQUIVALENCE (SYM(  1),SYN  1(1))
      EQUIVALENCE (SYM(  9),SYN  9(1))
      EQUIVALENCE (SYM( 17),SYN 17(1))
      EQUIVALENCE (SYM( 25),SYN 25(1))
      EQUIVALENCE (SYM( 33),SYN 33(1))
      EQUIVALENCE (SYM( 41),SYN 41(1))
      EQUIVALENCE (SYM( 49),SYN 49(1))
      EQUIVALENCE (SYM( 57),SYN 57(1))
      EQUIVALENCE (SYM( 65),SYN 65(1))
      EQUIVALENCE (SYM( 73),SYN 73(1))
      EQUIVALENCE (SYM( 81),SYN 81(1))
      EQUIVALENCE (SYM( 89),SYN 89(1))
      EQUIVALENCE (SYM( 97),SYN 97(1))
      EQUIVALENCE (SYM(105),SYN105(1))
      EQUIVALENCE (SYM(113),SYN113(1))
      EQUIVALENCE (SYM(121),SYN121(1))
      EQUIVALENCE (STB(  1),STC  1(1))
      EQUIVALENCE (STB(  9),STC  9(1))
      EQUIVALENCE (STB( 17),STC 17(1))
      EQUIVALENCE (STB( 25),STC 25(1))
      EQUIVALENCE (STB( 33),STC 33(1))
      EQUIVALENCE (STB( 41),STC 41(1))
      EQUIVALENCE (STB( 49),STC 49(1))
      EQUIVALENCE (STB( 57),STC 57(1))
      EQUIVALENCE (STB( 65),STC 65(1))
      EQUIVALENCE (STB( 73),STC 73(1))
      EQUIVALENCE (STB( 81),STC 81(1))
      EQUIVALENCE (STB( 89),STC 89(1))
      EQUIVALENCE (STB( 97),STC 97(1))
      EQUIVALENCE (STB(105),STC105(1))
      EQUIVALENCE (STB(113),STC113(1))
      EQUIVALENCE (STB(121),STC121(1))
      EQUIVALENCE (STB(129),STC129(1))
      EQUIVALENCE (STB(137),STC137(1))
      EQUIVALENCE (STB(145),STC145(1))
      EQUIVALENCE (STB(153),STC153(1))
      EQUIVALENCE (STB(161),STC161(1))
      EQUIVALENCE (STB(169),STC169(1))
      EQUIVALENCE (STB(177),STC177(1))
      EQUIVALENCE (STB(185),STC185(1))
      EQUIVALENCE (STB(193),STC193(1))
      EQUIVALENCE (STB(201),STC201(1))
      EQUIVALENCE (STB(209),STC209(1))
      EQUIVALENCE (STB(217),STC217(1))
      EQUIVALENCE (STB(225),STC225(1))
      EQUIVALENCE (STB(233),STC233(1))
      EQUIVALENCE (STB(241),STC241(1))
      EQUIVALENCE (STB(249),STC249(1))
      EQUIVALENCE (STB(257),STC257(1))
      EQUIVALENCE (STB(265),STC265(1))
      EQUIVALENCE (STB(273),STC273(1))
      EQUIVALENCE (STB(281),STC281(1))
      EQUIVALENCE (STB(289),STC289(1))
      EQUIVALENCE (STB(297),STC297(1))
      EQUIVALENCE (STB(305),STC305(1))
      EQUIVALENCE (STB(313),STC313(1))
      EQUIVALENCE (STB(321),STC321(1))
      EQUIVALENCE (STB(329),STC329(1))
      EQUIVALENCE (STB(337),STC337(1))
      EQUIVALENCE (STB(345),STC345(1))
      EQUIVALENCE (STB(353),STC353(1))
      EQUIVALENCE (STB(361),STC361(1))
      EQUIVALENCE (STB(369),STC369(1))
      EQUIVALENCE (STB(377),STC377(1))
      EQUIVALENCE (STB(385),STC385(1))
      EQUIVALENCE (STB(393),STC393(1))
      EQUIVALENCE (STB(401),STC401(1))
      EQUIVALENCE (STB(409),STC409(1))
      DATA SYN	1/     0,     8,    20,    26
     1,    31,	  36,	 43,	50/
      DATA SYN	9/    55,    62,    68,    82
     1,    91,	  98,	101,   104/
      DATA SYN 17/   111,   118,   125,   125
     1,   125,	 125,	125,   125/
      DATA SYN 25/   125,   125,   125,   125
     1,   125,	 125,	125,   125/
      DATA SYN 33/   125,   126,   136,   141
     1,   150,	 161,	174,   186/
      DATA SYN 41/   190,   195,   200,   207
     1,   212,	 215,	218,   221/
      DATA SYN 49/   223,   234,   240,   249
     1,   260,	 266,	276,   287/
      DATA SYN 57/   293,   311,   321,   326
     1,   332,	 336,	341,   345/
      DATA SYN 65/   354,   366,   373,   385
     1,   394,	 402,	409,   414/
      DATA SYN 73/   424,   430,   437,   443
     1,   448,	 452,	457,   463/
      DATA SYN 81/   473,   480,   492,   501
     1,   514,	 519,	526,   530/
      DATA SYN 89/   536,   544,   552,   559
     1,   564,	 567,	571,   577/
      DATA SYN 97/   583,   586,   597,   607
     1,   614,	 625,	635,   641/
      DATA SYN105/   654,   660,   668,   675
     1,   681,	 687,	697,   704/
      DATA SYN113/   714,   724,   735,   740
     1,   748,	 754,	762,   766/
      DATA SYN121/   772,   776,   785,   790
     1,   798,	 803,	811,   816/
      DATA STC	1/ 10842,  2088, 11276, -9686
     1, 10842,	8233,  2320,  5131/
      DATA STC	9/ 11044, -9686, 10842,  3080
     1, -9686,	2666,  7256, 18650/
      DATA STC 17/ 26668, -9716, 10842,  2584
     1, 10780, 23258, 10776,  6684/
      DATA STC 25/ -9718, 11336,  3112, 26842
     1,  2092, 22796, -9701, 10330/
      DATA STC 33/  6764, -9718, 11354,  8547
     1, 24872,	2065,  4945, 21260/
      DATA STC 41/ -9693, 11336,  3176, 10826
     1,  7256, 23258, 11304,  3080/
      DATA STC 49/ -9702, 10826, 22746, -9700
     1,  7256, 11368,  3144, 24794/
      DATA STC 57/ 20516, 18708, -9685, 10842
     1,  9057,	4945,-14630,   577/
      DATA STC 65/  2314, 20993, 12849,-14830
     1, 12641,	9075, 16838, 29489/
      DATA STC 73/ 24579, 21540,-14832,  2888
     1,  6932,	8217, 11305,   626/
      DATA STC 81/ 18630, 26924, 12328, 10545
     1,   836,	3083,-14844,   596/
      DATA STC 89/  2049,  8976, 12843,  6441
     1,-14844, 10849,-14798,  8563/
      DATA STC 97/   785, 29126,  4899,-14847
     1, 10826,	4196,  5216, 19142/
      DATA STC105/ 22570,-14820,  4361, 22726
     1,-14820,	 321, 13510, 18630/
      DATA STC113/ 12584, 11315,   780,  2049
     1,-14804,	 833, 12866,-14807/
      DATA STC121/ 12648, 11315,  2084,  1024
     1, 28870, 11316,  6938,  3092/
      DATA STC129/   259,-14840, 13123,  4120
     1,-14828,	 328,  3075,  8988/
      DATA STC137/ 12320,-14796,  7000,  3092
     1,   259,	8200, 13362, 28870/
      DATA STC145/ 11316,   273, 16838,  3075
     1,  6932,	8217, 12584, 11315/
      DATA STC153/  6948,  4185,   264,   710
     1, 11284, 12595,  8232,  7193/
      DATA STC161/ 16838, 20737,-14831,  2625
     1, 25106,-14814,  6211,-14797/
      DATA STC169/  5200,  8292, 16838, 12572
     1, 17094, 20994, 11300, 12595/
      DATA STC177/-14808,   323, 10248, 13105
     1,  5164,	6674,-14804, 12832/
      DATA STC185/  1060,  5200, 12486, 11315
     1,  6948, 23320,  3092,	 3/
      DATA STC193/ 27846, 12595,  2088,   769
     1,-14836,	3075, 13100, 28976/
      DATA STC201/-14847, 13360,  6235,  1088
     1, 12486, 23092,-14824, 13164/
      DATA STC209/ 10289,   264,  5124,-14829
     1, 22576, 29724,-14844,   833/
      DATA STC217/ 12866, 13169, 18630,   769
     1, 13324, 12486,  6260,-14844/
      DATA STC225/ 16432,-14844,  6704,  1076
     1, 12486,	3176,  1140, 18630/
      DATA STC233/   769, 11276, 12595,  2088
     1, 12486, 11315,  6948,-14824/
      DATA STC241/   328,  3075, 13100, 10289
     1, 21000,-14844, 13104,  9260/
      DATA STC249/  6171,  1114, 18630,   769
     1,  5132,	6427, 10272, 13105/
      DATA STC257/-14804, 13424,   626, 28870
     1,   264,	3075,-14796,   624/
      DATA STC265/-14796,   368,   794,-14796
     1, 11272, 28724,  3112,-14844/
      DATA STC273/ 10352,   538, 11380,-14822
     1, 13424,	2092,  1024, 29382/
      DATA STC281/    48,-15358, 12356,   710
     1, 12338, 22724,  7210,  2666/
      DATA STC289/ 19142, 10776,  7256, 28870
     1,-14821,	6216,  8737,   532/
      DATA STC297/  2049,  1124, 12486,  8784
     1,  7203,	 780,  4098, 25798/
      DATA STC305/  6177,   264,-14844,  1140
     1,  8788,	6177,	264,  5122/
      DATA STC313/ 20678,  7188,  8483,  2072
     1,   769, 17094, 13098,  7001/
      DATA STC321/ 16838,  3075, 25652, 12594
     1,  6184,	4625,-14812, 24624/
      DATA STC329/  7203,-14844,   833,  6722
     1, 27161,-14806,	577,  8971/
      DATA STC337/ 13171, 12486,  2156,  1114
     1, 16838, 16899, 12594,  8390/
      DATA STC345/  8536,   538,  9050,  1052
     1,  8390,	8536,  7203,-14844/
      DATA STC353/  6216,  8993,  3100,   259
     1,-14840, 24624, 13106,  7212/
      DATA STC361/  4627,-14816, 13380, 12900
     1, 10289,	4376,  9234,  8390/
      DATA STC369/  8784,-14812,  3075,  4371
     1,  8472,-14812,  9057,  2674/
      DATA STC377/-14845,  2144,   769, 25612
     1,-14844,	 608,-14812,   352/
      DATA STC385/   786,-14812, 24612,-14844
     1,   833, 13324,  6256,  5137/
      DATA STC393/ 24774,    36,-14844,  2370
     1,  6161, 10529,-15054,  4674/
      DATA STC401/ 12898, 16838,  4618,  8731
     1, 12586, 20677,  6937,-14812/
      DATA STC409/ 13360,     4, 28724,-14844
     1,     0,	   0,	  0,	 0/
C
	DIMENSION ROT(4,4),TEMP(4,4),V(4)
	DIMENSION XX(2),YY(2),ZZ(2)
C
	IOPT=IOP		! VISIBILITY FLAG
	XA1=XA
	YA1=YA
	ZA1=ZA
	IF (XA.EQ.999..OR.YA.EQ.999..OR.ZA.EQ.999.) THEN !CONTINUE STRING OUTPUT
		XA1=X
		YA1=Y
		ZA1=Z
	ELSE
C
C	SET UP 3-D CONVERSION MATRIXES
C
		CALL ROTEM(1,PG,ROT)
		CALL ROTEM(3,PB,TEMP)
		CALL MATMUL4(ROT,ROT,TEMP)
		CALL ROTEM(2,PA,TEMP)
		CALL MATMUL4(ROT,ROT,TEMP)
	ENDIF
	X=XA1
	Y=YA1
	Z=ZA1
C
	H=HI/7.0
	IF (NC.NE.-2) GOTO 4
	IP=2
	IF (IOPT.EQ.0) CALL PLT3DH(XA1,YA1,Z,IP)
	XX(1)=XA1
	YY(1)=YA1
	ZZ(1)=Z
4	IP=3
	IF (IOPT.EQ.0) CALL PLT3DH(X,Y,Z,IP)
	XX(1)=X
	YY(1)=Y
	ZZ(1)=Z
	IF (NC) 10,20,30
10	NC0=0
	JCH=ASC(1)+0.1
	GOTO 50
20	RETURN
30	NC0=NC
	I0=0
40	NC0=NC0-1
	IF (NC0.LT.0) GOTO 20
	I0=1+I0
C
C	PXPCGT GET THE ASCII VALUE OF THE I0TH CHARACTER IN ASC
C	AND PUTS IT INTO JCH
C
	CALL PXPCGT(ASC,I0,JCH)
50	IF (JCH.LT.0.OR.JCH.GE.128) JCH=32
	I3=SYM(1+JCH)
	IF (JCH.GE.32) GOTO 55
C	CENTERED PLOTTING SYMBOLS
	IP=3
	I1=0
	V(1)=-H*2.
	V(2)=-H*3.
	V(3)=0.
	CALL MTV4(V,ROT,V)
	X=V(1)+X
	Y=V(2)+Y
	Z=V(3)+Z
	XW1=0.
	YW1=0.
	XW2=0.
	YW2=0.
	GOTO 110
55	IF (JCH.NE.112) GOTO 60
	XW2=0.
	YW2=-H*2.
	V(1)=XW2
	V(2)=YW2
	V(3)=0.0
	CALL MTV4(V,ROT,V)
	IF (IOPT.EQ.0) CALL PLT3DH(V(1)+X,V(2)+Y,V(3)+Z,3)
	XX(1)=V(1)+X
	YY(1)=V(2)+Y
	ZZ(1)=V(3)+Z
60	IP=2
	I3=1+I3
	I4=(I3-1)/2+1
	IV=STB(I4)
	IJ=MOD(I3-1,2)*8
	I1=IBITS(IV,IJ,8)
	I2=MOD(I1,64)
	IF (I1.GE.64) IP=3
	IY=I2/8
	IX=I2-8*IY
	X1=IX
	Y1=IY
	XW1=H*X1
	YW1=H*Y1
	XW2=XW1
	YW2=YW1
	IF ((JCH.EQ.103.OR.JCH.EQ.112.OR.JCH.EQ.113.OR.
     $		JCH.EQ.121).AND.(I1.LT.192)) THEN
		XW2=XW2
		YW2=YW2-H*2.
	ENDIF
110	CONTINUE
	V(1)=XW2
	V(2)=YW2
	V(3)=0.0
	CALL MTV4(V,ROT,V)
	IF (IOPT.EQ.0) THEN
		CALL PLT3DH(V(1)+X,V(2)+Y,V(3)+Z,IP)
	ELSE
		IF (IP.EQ.3) THEN
			XX(1)=V(1)+X
			YY(1)=V(2)+Y
			ZZ(1)=V(3)+Z
		ELSE
			XX(2)=V(1)+X
			YY(2)=V(2)+Y
			ZZ(2)=V(3)+Z
			IER=0
			NN=2
			IF (ABS(XX(1)-XX(2)).GT.0.01.OR.ABS(YY(1)-YY(2))
     $			  .GT.0.01.OR.ABS(ZZ(1)-ZZ(2)).GT.0.01)
     $     			CALL SKETCH(XX,YY,ZZ,NN,IER)
			IF (IER.NE.0) THEN
	CALL CTERM(1)
	WRITE(*,1400)
1400	FORMAT(' *** WARNING: SYM3DH SKETCH ERROR')
	CALL CTERM(-1)
			ENDIF
			XX(1)=XX(2)
			YY(1)=YY(2)
			ZZ(1)=ZZ(2)
		ENDIF
	ENDIF
	IF (I1.LT.192) GOTO 60
	V(1)=XW1
	V(2)=YW1
	V(3)=0.0
	CALL MTV4(V,ROT,V)
	X=V(1)+X
	Y=V(2)+Y
	Z=V(3)+Z
	GOTO 40
	END
C
C
	SUBROUTINE NUM3DH(XA,YA,ZA,PA,PB,PG,HGHT,Z,F0,IOP)
C
C	CREATED: DGL 1-OCT-1984
C	MODIFIED: DGL FEB 1986 -- ADDED IOP WITH SKETCH CALL
C	
C	PLOTS A FLOATING POINT NUMBER IN 3-D USING SYM3DH
C
C	XA,YA,ZA COORDINATES OF STRING LOWER-LEFT CORNER
C		  (999.,999.,999.) CONTINUE FROM LAST POINT
C	PA,PB	SPECIFIES THE PROJECTED ANGLES (IN DEGS) RELATIVE TO THE
C		 X-Y PLANE AND THE X-Z PLANE OF THE LINE ALONG THE
C		 BASE OF THE PLOTTED NUMBER SEQUENCE
C	PG	SPECIFIES THE ROTATION OF THE PLOTTED SYMBOLS AROUND
C		THE LINE SPECIFIED BY PA,PB
C	HGHT	HEIGHT OF THE PLOTTED NUMBER
C	Z	FLOATING POINT NUMBER TO BE PLOTTED
C	F0	PLOTTING FORMAT (Fn.j) [DEFAULT 0.0]
C		  n = TOTAL NUMBER OF SPACES TO USE (INCLUDING SIGN AND D.P.)
C			[MAX 18 CHARACTERS WIDE]
C		  j = DIGITS TO RIGHT OF DECIMAL POINT (2 DIGITS)
C		  F0 < 0 SPECIFIES EXPONENTIAL NOTATION
C		  F0 = PLOT INTEGER PORTION ONLY WITH NO D.P. (FREE FORMAT)
C		  F0 > 1000 PLOT INTEGER PORTION IN FIXED FORMAT
C		  IF n=0 THEN J IS NUMBER OF DIGITS TO RIGHT OF DECIMAL POINT
C		  WHEN Z OVERFLOWS THIS FORMAT, SPACE IS FILLED WITH ASTERICKS
C	IOP	FLAG TO INDICATE SKETCH OR PLT3DH USAGE
C		IOP=0: PLT3DH
C		IOP=1: SKETCH
C NOTE:
C	F0 IS SET UP TO HANDLE AN INTEGER ARGUMENT (MACHINE DEPENDENT)
C		(NUMBER OF DIGITS TO RIGHT OF DECIMAL POINT)
C
	BYTE B(18)
	EQUIVALENCE (NN,RR)
	ND=-1
	XT=XA
	YT=YA
	ZT=ZA
	PAT=PA
	PBT=PB
	PGT=PG
	HG=HGHT
	FA=F0
	IOPT=IOP
	IF (FA.NE.0.0) GOTO 3
	RR=FA
	IF (NN.NE.0) GOTO 2		! CHECK FOR INTEGER F0 INPUT
    	ND=0
	FA=1.0
	GOTO 5
  1	NN=AMOD(FA,1000.)		! PLOT FORMATED INTEGER 
	FA=0.
	GOTO 5
  2	FA=1.0
	IF (NN.LT.0) FA=0.0
	ND=NN
	NN=0
	GOTO 5
  3	IF (ABS(FA).GT.1.E30)  THEN
		RR=FA
		GOTO 2			! INTEGER -1 INPUT
	ENDIF
	IF (FA.GT.999) GOTO 1
	F=ABS(FA)*1.000002		! GET PLOT FORMAT
	NN=F
	F=(F-NN)*100.
	ND=F
  5	IF (ND.GT.17) ND=ND/10		! SOLVE SIMPLE ERRORS
	IF (NN.EQ.0) THEN		! DIGITS TO LEFT OF DECIMAL POINT
		NN=ND+2
		IF (Z.EQ.0.AND.FA.EQ.0.0) NN=1
		IF (Z.NE.0.) THEN
			ALG=ALOG10(ABS(Z))
			IF (ALG.LT.0.0) ALG=0.
			NN=ND+2+ALG
		ENDIF
		IF (Z.LT.0.0) NN=NN+1
		IF (FA.LT.0.0) NN=NN+4
	ENDIF
	IF (ND.GT.NN) GOTO  60		! FORMAT ERROR
	IF (NN.GT.18) NN=18		! MAX CHARACTERS
	IF (FA) 20,30,40
 20	ENCODE(18,21,B,ERR=90) Z
 21	FORMAT(E<NN>.<ND>)
	GOTO 50
 30	I=Z
	ENCODE(18,31,B,ERR=90) I
 31	FORMAT(I<NN>)
	GOTO 50
 40	ENCODE(18,41,B,ERR=90) Z
 41	FORMAT(F<NN>.<ND>)
 50	CALL SYM3DH(XT,YT,ZT,PAT,PBT,PGT,HG,B,NN,IOPT)
 60	RETURN
 90	DO 95 I=1,18
		B(I)='*'
		IF (I.EQ.NN-ND) B(I)='.'
 95	CONTINUE
	GOTO 50
	END
C
C *************************************************************************
C
C	ROUTINES FOR PLOTTING IN 3-D WITH HIDDEN LINE REMOVAL
C	(MODIFIED FROM ORIGINALS FROM COSMIC)
C	(LAST REVISED: FEB. 15, 1986  DGL)
C
C	TO USE HIDDEN LINE REMOVAL:
C		1. SET UP WORKING ARRAY IN COMMON IN YOUR MAIN PROGRAM:
C			COMMON /GO/ISIZE,WORK(ISIZE)
C		     WHERE ISIZE >= (25 + 5*MNE + 4*(MNE+2+2*MNH))*NPOLYS
C		     WHERE NME=MAX NUMBER OF EDEGES PER POLY,
C			   MNH=MAX NUMBER OF HOLES PER POLY,
C			   NPOLYS=MAX NUMBER OF SCENE ELEMENTS
C		1. CALL SKETCH WITH ALL BUT LAST POLYGON USING NC=0
C		2. CALL SKETECH WITH LAST POLYGON USING NC=1
C
C ***********************************************************************
C
	SUBROUTINE SKETCH(X,Y,Z,NPP,NC)
C
C     THIS ROUTINE SETS UP PEN MOTION INDICATORS.
C     FOR THE HIDDEN LINE REMOVAL ALGORITHM
C
C	X,Y,Z	(R)	ARRAYS X(NP) CONTAINING POLYGON EDGES
C	NPP	(I)	NUMBER OF POINTS ON POLYGON
C	NC	(I)	FLAG.  SHOULD BE:
C				SET TO 0 FOR ALL BUT THE LAST CALL
C				SET TO 1 FOR LAST CALL
C			NOTE: SETTING NC=FLAG+10 WILL CAUSE ONLY THE FIRST
C			      VALUE OF Z TO BE REPLICATED FOR ALL X,Y PAIRS
C			RETURNED VALUE:
C				=  0 NORMAL RETURN
C				= -1 ICORRECT STORAGE ALLOCATION IN /GO/ 
C				= -2 DISTANCE FROM VIEW INCORRECT (DV MODFIED)
C				= -3 DV AND /GO/ ERRORS
C				= -4 TO MANY EDGES FOR ONE POLYGON
C
	DIMENSION X(1),Y(1),Z(1)
C
	LOGICAL ZFILL
	DIMENSION X1(160),Y1(160),Z1(160) ! 160 = MAXIMUM NUMBER OF POLY SIDES
	COMMON /CHIDE/AL(3),V0(3),VSF,SCX,YAW,ROL,PIT,LZ,VP,JJJ
C
	IF (NC.LT.0) RETURN	! DO NOT PROCESS IF PRIOR ERROR HAS OCCURED
	IF (NPP.GT.LZ) THEN	! TO MANY EDGES ON POLYGON
		IERR=-4
		RETURN
	ENDIF
	ZFILL=.TRUE.
	IF (NC/10.EQ.1) THEN	! SPECIAL Z VALUE REPLICATION
		ZFILL=.FALSE.
		NC=MOD(NC,10)
	ENDIF
	NP=NPP
	L=NP
	LI=NP
	IF (L.LE.2) GOTO 50
	LX=1
	NPX=NP
    1	NPX=NPX-1
	I=LX
	DO 8 M=I,NPX
		RX=0
		A=X(M+1)-X(M)
		B=Y(M+1)-Y(M)
		IF (ZFILL) THEN
			C=Z(M+1)-Z(M)
		ELSE
			C=0.
		ENDIF
		IF (A.NE.0.) GOTO 8
		IF (B.NE.0.) GOTO 8
		IF (C.NE.0.) GOTO 8
		IX=M
		IX1=NPX
		DO 4 MX=IX,IX1
			X(MX)=X(MX+1)
			Y(MX)=Y(MX+1)
			IF (ZFILL) Z(MX)=Z(MX+1)
    4	CONTINUE
	RX=1
	LX=M
	IF (LX.EQ.NPX) GOTO 10
	GOTO 1
    8	CONTINUE
   10	CONTINUE
	IF (RX.EQ.1.) NPX=NPX-1
	NP=NPX+1
	LI=NP
	IF (NP.LE.2) GOTO 50
	IX=0
	M1=0
	M=1
	IS=NP-1
20	CONTINUE
	M=M+IX
	M1=M1+IX+1
	IF (M-1.EQ.LI) GOTO 70
C
C     SEARCH FOR MATCHING COORDINATES.
C
	DO 40 J=M,IS
		T=X(J+1)-X(M)
		U=Z(J+1)-Z(M)
		IF (ZFILL) THEN
			W=Y(J+1)-Y(M)
		ELSE
			W=0.
		ENDIF
		IF (T.NE.0.) GOTO 40
		IF (W.NE.0.) GOTO 40
		IF (U.NE.0.) GOTO 40
		NP=NP+1
	
C     MATCH FOUND.....STORE COORDINATES AND SET SWITCH TO LIFT PEN
C     AND/OR END SET.
C
		IX=J+2-M
		IX1=J-IS+1
		DO 30 IK=1,IX
			X1(M1-1+IK)=X(M-1+IK)
			Y1(M1-1+IK)=Y(M-1+IK)
			IF (ZFILL) THEN
				Z1(M1-1+IK)=Z(M-1+IK)
			ELSE
				Z1(M1-1+IK)=Z(1)
			ENDIF
   30		CONTINUE
		Z1(M1+IX)=-ISIGN(1,IX1)*9999.
		GOTO 20
   40	CONTINUE
   50	CONTINUE
	DO 60 J=1,LI
		X1(J)=X(J)
		Y1(J)=Y(J)
		IF (ZFILL) THEN
			Z1(J)=Z(J)
		ELSE
			Z1(J)=Z(1)
		ENDIF
   60	CONTINUE
	NP=NP+1
	Z1(NP)=-9999.
   70	CONTINUE
	IF (NP.GT.LZ) THEN	! TOO MANY EDGES ON ONE POLYGON
		IERR=-4
		RETURN
	ENDIF
	CALL LINHID(X1,Y1,Z1,NP,NC)
	NP=L
C
	RETURN
	END
C
      SUBROUTINE LCOEF(X,Y,Z,XXX,JXX,NC,NS,CCC,LZ)
C
C     THIS ROUTINE DETERMINES EQUATION OF LINES AND PLANES.
C
      DIMENSION CCC(1),XXX(1),X(1),Y(1),Z(1)
      DIMENSION COE(8)
      DIMENSION IB(5)
      DIMENSION ZZ(150),ZZZ(150)
      COMMON/GO3/L0,L1,L00,L01,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13
C
      LE=0
      JA=L13+(JXX-1)*LZ
      JF=L12+(JXX-1)*5
      I=0
      J=1
   10 CONTINUE
C
C     SEARCH FOR MATCHING COORDINATES.
C
      I=I+1
      T=X(I+1)-X(I)
      S=Y(I+1)-Y(I)
      U=Z(I+1)-Z(I)
      IF(T.NE.0.)GOTO 20
      IF(S.NE.0.)GOTO 20
      IF(U.NE.0.)GOTO 20
C
C     MATCH FOUND.....PROCEED IF LIST IS NOT EXHAUSTED.
C
      I=I+2
   20 CONTINUE
      IF(I.GT.NS)GOTO 70
C
C     DETERMINE EQUATION OF LINE-SEGMENTS.
C
      T=X(I+1)-X(I)
      T1=Y(I+1)-Y(I)
      IF((T1.EQ.0.).AND.(T.EQ.0.))GOTO 10
      IF(T.NE.0.)GOTO 30
   29 CONTINUE
      CCC(J+JA)=0
      CCC(J+1+JA)=1
      CCC(J+2+JA)=-X(I)
      GOTO 40
   30 CONTINUE
      CCC(J+JA)=1
      E=(Y(I+1)-Y(I))/(X(I+1)-X(I))
      F=(E*X(I))-Y(I)
      IF(ABS(E).LT.100.)GOTO 220
      YO=E*(X(I+1))-F
      G=ABS(Y(I+1)-YO)
      IF(G.GT..00001)GOTO 29
  220 CONTINUE
      CCC(J+1+JA)=-E
      CCC(J+2+JA)=F
   40 CONTINUE
      IF(CCC(J+JA).NE.0.)GOTO 50
      CCC(J+3+JA)=Y(I)
      CCC(J+4+JA)=Y(I+1)
      GOTO 60
   50 CONTINUE
      CCC(J+3+JA)=X(I)
      CCC(J+4+JA)=X(I+1)
   60 CONTINUE
      J=J+5
      LE=LE+1
      ZZ(LE)=Z(I)
      ZZZ(LE)=Z(I+1)
      IF(LE.GT.3)GOTO 10
      IB(LE)=I
      GOTO 10
   70 CONTINUE
C
C     DETERMINE EQUATION OF PLANE.
C
      J=(J-1)/5
      XXX(JF+5)=J
      IF(NS.LE.3)GOTO 120
      K1=1
      K2=2
      K3=3
      A1=X(K3)-X(K1)
      B1=Y(K3)-Y(K1)
      C1=Z(K3)-Z(K1)
      A2=X(K2)-X(K1)
      B2=Y(K2)-Y(K1)
      C2=Z(K2)-Z(K1)
      COE(1)=B1*C2-B2*C1
      COE(2)=C1*A2-C2*A1
      COE(3)=A1*B2-A2*B1
      COE(4)=COE(1)*X(1)+COE(2)*Y(1)+COE(3)*Z(1)
      COE(4)=-COE(4)
      DO 110 J=1,4
  110 XXX(JF+J)=COE(J)
      IF(COE(3).NE.0.) GOTO 140
      J=1
      DO 25 K=1, LE
      CCC(JA+J)=ZZ(K)
      CCC(JA+J+1)=ZZZ(K)
      J=J+5
   25 CONTINUE
      IF(COE(1).NE.0.)I=1
      IF(COE(2).NE.0.)I=2
      P=COE(I)
      DO 26 K=1,4
   26 XXX(JF+K)=XXX(JF+K)/P
      GOTO 140
  120 CONTINUE
      XXX(JF+5)=1
      DO 130 IX=1,2
  130 XXX(JF+IX)=Z(IX)
      XXX(JF+3)=0
  140 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE LINHID(X,Y,Z,NP,NC)
C
C     THIS ROUTINE IS THE EXECUTIVE ROUTINE FOR HIDDEN LINE REMOVAL
C
      INTEGER XIND(1),RCT(100)
      DIMENSION RRX(20)
      DIMENSION XXX(1),CCC(1)
      DIMENSION TGM(1),IN(1),ZM(1),ZMI(1),
     1TGMT(1),TGI(1),NNO(1),RV(1),RVI(1),NOCT(1)
      DIMENSION NGX(15)
      DIMENSION X21(200),Y21(200),Z21(200),IIA(200)
      DIMENSION XI(1000,1),YI(1000,1),ZI(1000,1),DI(1000)
      DIMENSION U(6),V(6),W(6),X(1),Y(1),Z(1),X1(10),Y1(10),Z1(10),H(15)
      DIMENSION ZMIN(1),YMIN(1),IN1(1),IN2(1),ICOUNT(150)
      DIMENSION XE(150),YE(150),XU(150),YU(150)
      DIMENSION I2(2),I3(2)
      DIMENSION COORD(1)
      DIMENSION SNDT(1)
      DIMENSION IND(1)
      DIMENSION REX(20)
      DIMENSION NEH(1),KEEP(1)
      DIMENSION IBEG(100),IEND(100),ICT(100),ICCT(100)
C
      COMMON /GO/ICORE,WORK(1)
      COMMON /GO3/L0,L1,L00,L01,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13
      COMMON /CHIDE/AL(3),V0(3),VSF,SCX,YAW,ROL,PIT,LZ,VP,JJJ
      COMMON /HEDG/JAT,ME,JT,VX,VX1,VX2,VX3,NN
C 
      EQUIVALENCE(WORK(1),XXX(1),CCC(1),IN(1),TGMT(1))
      EQUIVALENCE(WORK(1),ZM(1),ZMI(1),NNO(1),
     1	TGM(1),TGI(1),RV(1),RVI(1),NOCT(1))
      EQUIVALENCE(WORK(1),IN1(1),IN2(1),YMIN(1),ZMIN(1))
      EQUIVALENCE(WORK(1),COORD(1))
      EQUIVALENCE(WORK(1),SNDT(1))
      EQUIVALENCE (WORK(1),IND(1),XIND(1))
      EQUIVALENCE(WORK(1),KEEP(1),NEH(1))
C 
      IF(VP.LT.0.)GOTO 20
      HXX=.015
      AVA=0
      HX1=.001
      LC=10**6
      IXXX=0
      IF(SCX.LT.0.)IXXX=1
      SCX=ABS(SCX)
      SCT=1
      VP=VP/SCT
C
C     INITIALIZE VARIABLES.
C
      VPX=ABS(VP)
      LZ=LZ*5		! TEMPORARILY UPSCALE
      SW1=0
      SW=0
      IDAV=0
C
C     CALCULATE MAXIMUM ALLOWABLE ELEMENTS.
C
      IABC=ICORE/(25+LZ+4*JJJ)
      NCC=IABC
      L5=0
      L6=NCC
      L7=2*NCC
      L8=3*NCC
      L2=4*NCC
      L3=5*NCC
      L4=6*NCC
      L00=7*NCC
      L01=8*NCC
      L1=9*NCC
      L0=10*NCC
      L9=11*NCC
      L10=12*NCC
      L11=13*NCC
      L15=14*NCC
      L16=15*NCC
      L17=16*NCC
      L18=19*NCC
      L12=20*NCC
      L13=25*NCC
      L14=L13+LZ*NCC
      DO 10 J=1,NCC
      RVI(L8+J)=10**6
      TGM(L5+J)=10**6
      RV(L7+J)=-RVI(L8+J)
      TGI(L6+J)=-TGM(L5+J)
      NOCT(L9+J)=0
      ZM(L2+J)=RV(L7+J)
      ZMI(L3+J)=RVI(L8+J)
      XIND(L16+J)=0
      IND(L15+J)=J
      KEEP(L18+J)=0
   10 CONTINUE
      IK=0
      IKT=0
      KR=JJJ
      PI=3.1416/180.
      U(6)=SCX
      V(6)=SCX
      VP=-VP
C
C     STORE EULERIAN ANGLES.
C
      XX=YAW*PI
      YY=ROL*PI
      ZZ=PIT*PI
      COSY=COS(YY)
      SINY=SIN(YY)
      COSZ=COS(ZZ)
      SINZ=SIN(ZZ)
      COSX=COS(XX)
      SINX=SIN(XX)
   20 CONTINUE
      NT=NP-1
       IKK=IK+1
       IK=IK+1
C
C     SET ERROR CODES, IF NECESSARY.
C
      IF (IKK.LE.IABC) GOTO 30
      SW=1				! MEMORY ERROR
   30 CONTINUE
      IF (NC.EQ.0) GOTO 40
      IDAV=1				! LAST PASS FLAG
      NC=-SW1
      IF (SW.EQ.0.) GOTO 50
      NC=-(SW+SW1)
      ICORES=(25+LZ+4*JJJ)*IKK		! COMPUTE STORAGE SPACE SIZE
      IF (ICORES.GT.ICORE) RETURN	! DO NOT RESET--STOP ON ERROR
   40 CONTINUE
   50 CONTINUE
      DO 60 J=1,NP
	X21(J)=X(J)
	Y21(J)=Y(J)
	Z21(J)=Z(J)
   60 CONTINUE
C
C     STORE COORDINATES AND SET PEN POSITION WHENEVER ABS(Z)=9999.
C
      DO 70 J=1,NT
      IIA(J)=0
      IF(Z21(J).NE.9999.)GOTO 70
      IIA(J)=1
      IXU=J-2
      IBB=J-ISIGN(1,IXU)
      X21(J)=X21(IBB)
      Y21(J)=Y21(IBB)
      Z21(J)=Z21(IBB)
   70 CONTINUE
      IIA(NP)=1
      Z21(NP)=Z21(NT)
      Y21(NP)=Y21(NT)
      X21(NP)=X21(NT)
      JXX=IKK
      I=1
      VL=ABS(VP)
C
C     LOOP THAT DOES THE THREE DIMENSIONAL TRANSFORMATION ON THE
C     COORDINATES.
C
      JV=L14+(IKK-1)*4*JJJ
      JT=1
      DO 90 J=1,NP
      XJ=VSF*X21(J)/SCT+V0(1)
      YJ=VSF*Y21(J)/SCT+V0(2)
      ZJ=VSF*Z21(J)/SCT+V0(3)
      U(I)=ZJ*(COSY*SINX)+XJ*(COSY*COSX)-YJ*SINY
      TW=YJ*COSY*COSZ
      TZ=XJ*(SINZ*SINX+SINY*COSZ*COSX)
      TY=ZJ*(-SINZ*COSX+SINY*COSZ*SINX)
      V(I)=TZ+TW+TY
      PT=YJ*COSY*SINZ
      PK=ZJ*(COSZ*COSX+SINY*SINZ*SINX)
      PS=XJ*(-COSZ*SINX+SINY*SINZ*COSX)
      ZJ=PK+PS+PT
      IF(ZJ.LT.VL)GOTO 80
      SW1=2
      VPX=AMAX1(ZJ,VPX)
      VPX=VPX+.5
   80 CONTINUE
      T=SW+SW1
      IF(T.NE.0.)GOTO 90
C
C     CALCULATES PERSPECTIVE BASED ON VALUE VP(DV) FROM CALLING
C     PROGRAM.
C
      HH=VL/(VL-ZJ)
      X21(J)=U(I)*HH
      Y21(J)=V(I)*HH
      Z21(J)=ZJ*HH
C
C     CALCULATES MAX/MIN VALUES OF EACH ELEMENT ON THE X,Y,Z DIMENSION
C
      RV(L7+JXX)=AMAX1(RV(L7+JXX),Y21(J))
      RVI(L8+JXX)=AMIN1(RVI(L8+JXX),Y21(J))
      TGI(L6+JXX)=AMAX1(TGI(L6+JXX),X21(J))
      TGM(L5+JXX)=AMIN1(TGM(L5+JXX),X21(J))
      ZM(L2+JXX)=AMAX1(ZM(L2+JXX),Z21(J))
      ZMI(L3+JXX)=AMIN1(ZMI(L3+JXX),Z21(J))
      COORD(JT+JV)=X21(J)
      COORD(JT+JV+1)=Y21(J)
      COORD(JT+JV+2)=Z21(J)
      COORD(JT+3+JV)=IIA(J)
      JT=JT+4
   90 CONTINUE
      IF (IDAV.EQ.1) VP=VPX
      IF (T.NE.0.) GOTO 400
      NOCT(L9+IKK)=NOCT(L9+IKK)+NP
      NS=NP
      AVA=AVA+(TGI(L6+JXX)-TGM(L5+JXX))*(RV(L7+JXX)-RVI(L8+JXX))
      IF(IXXX.EQ.1)GOTO 95
C
C     CALL ROUTINE WHICH CALCULATES BOTH THE EQUATIONS OF THE LINE
C     SEGMENTS AND POLYGONS.
C
      CALL LCOEF(X21,Y21,Z21,XXX,JXX,NC,NS,CCC,LZ)
C
C     CHECKS TO SEE IF ALL ELEMENTS(SETS) HAVE BEEN PASSED.
C
   95 CONTINUE
      IF (IDAV.EQ.1) GOTO 100
      GOTO 400
  100 CONTINUE
      AVA=AVA/IKK
      DO 1301 J=1,100
      ICCT(J)=0
      ICT(J)=0
      RCT(J)=J-1
      IBEG(J)=1
      IEND(J)=0
 1301 CONTINUE
      IAUG=50+(IKK/10000)*2
      AMAXX=-9999.
      AMAXY=-9999.
      AMINX=9999.
      AMINY=9999.
      DO 1400 J=1,IKK
      AMAXX=AMAX1(AMAXX,TGI(L6+J))
      AMAXY=AMAX1(AMAXY,RV(L7+J))
      AMINX=AMIN1(AMINX,TGM(L5+J))
      AMINY=AMIN1(AMINY,RVI(L8+J))
 1400 CONTINUE
      TMAX=(AMAXX-AMINX)*(AMAXY-AMINY)
      IBL=TMAX/AVA
      IBL=IBL/4
C
C     DETERMINES THE NUMBER OF GRID POINTS IN THE GRID.
C
      EN=IKK
      K=(ALOG(EN)/ALOG(2.))+.01
      K=K+IAUG
      K=MIN0(K,IBL)
      IF(K.LE.1)K=1
      T=K
      R=(T**.5)
      KS=R+.5
      S=T/KS
      MS=S+.5
      N=KS*MS
      MND=N+1
      XMD=MND
      T=3./(MND-1)
      IGY=T*IKK
      K=KS
      K1=MS
      CRX=(AMAXX-AMINX)/K
      CRY=(AMAXY-AMINY)/K1
C
C     DETERMINES THE RELEVANT ELEMENTS VIA THE GRID BLOCKS.
C
      DO 3 J=1,IKK
      IA=0
      XMAT=TGI(L6+J)
      XMIT=TGM(L5+J)
      YMAT=RV(L7+J)
      YMIT=RVI(L8+J)
      M=0
      DO 1 I=1,K1
      DO 2 L=1,K
      M=M+1
      S=XMAT-((L-1)*CRX+AMINX)
      S1=XMAT-(L*CRX+AMINX)
      R=XMIT-((L-1)*CRX+AMINX)
      R1=XMIT-(L*CRX+AMINX)
      A=YMAT-((I-1)*CRY+AMINY)
      A1=YMAT-(I*CRY+AMINY)
      B=YMIT-((I-1)*CRY+AMINY)
      B1=YMIT-(I*CRY+AMINY)
      IF((S.LE.0.).OR.(R1.GE.0.))GOTO 2
      IF((A.LE.0.).OR.(B1.GE.0.))GOTO 2
      IF((S*S1.GT.0.).OR.(R*R1.GT.0.))GOTO 4
      IF((A*A1.GT.0.).OR.(B*B1.GT.0.))GOTO 4
      XIND(L16+J)=M
      GOTO 3
    4 CONTINUE
      IA=IA+1
      IF(IA.LE.4)GOTO 8000
      XIND(J+L16)=0
      GOTO 8001
 8000 CONTINUE
      XIND(L16+J)=XIND(L16+J)+M*(MND**(IA-1))
 8001 CONTINUE
      IF(ICCT(M).LT.0)GOTO 2
      ICCT(M)=ICCT(M)+1
      JK=(M-1)*IGY+ICCT(M)+L17
      NEH(JK)=J
      IF(ICCT(M).GE.IGY)ICCT(M)=-1
    2 CONTINUE
    1 CONTINUE
    3 CONTINUE
      CALL VSRT1(XIND(L16+1),IK,IND(L15+1))
      SW=0
      L=1
      DO 5 I=1,IKK
   11 CONTINUE
      IF(XIND(L16+I).NE.RCT(L))GOTO 6
      SW=SW+1
      IF(SW.EQ.1.)LT=I
      ICT(L)=ICT(L)+1
      GOTO 5
    6 CONTINUE
      IF(SW.NE.0.)GOTO 8
      L=L+1
      GOTO 11
    8 CONTINUE
      IBEG(L)=LT
      IEND(L)=LT+ICT(L)-1
      SW=0
      IF(XIND(L16+I).GE.MND)GOTO 2110
      L=L+1
      GOTO 11
    5 CONTINUE
      IBEG(L)=LT
      IEND(L)=LT+ICT(L)-1
 2110 CONTINUE
      DO 2111 J=1,IKK
      SNDT(L4+J)=IND(L15+J)
 2111 CONTINUE
      CALL VSRTR(SNDT(L4+1),IK,XIND(L16+1))
      EN=IKK
      IGX=(ALOG(EN)/ALOG(2.))+1.
      DO 105 J=1,IGX
      RRX(J)=2**(IGX-J)
  105 CONTINUE
      U(6)=SCX
      V(6)=SCX
      W(6)=SCX
      IKT=NC
      T=AMINY
      T1=AMINX
      V(5)=T
      U(5)=T1
      IJ=0
      X1(3)=U(5)
      Y1(3)=V(5)
      X1(4)=U(6)
      Y1(4)=V(6)
      X1(4)=X1(4)/SCT
      Y1(4)=Y1(4)/SCT
      DO 115 J=1,IKK
      IN(L11+J)=J
      IN1(L0+J)=J
      IN2(L00+J)=J
      TGMT(L10+J)=TGM(L5+J)
      YMIN(L1+J)=RVI(L8+J)
      ZMIN(L01+J)=ZM(L2+J)
  115 CONTINUE
C
C     CALL ROUTINE WHICH WILL SORT ON X,Y AND Z.
C
      CALL VSRTR(TGMT(L10+1),IK,IN(L11+1))
      CALL VSRTR(YMIN(L1+1),IK,IN1(L0+1))
      CALL VSRTR(ZMIN(L01+1),IK,IN2(L00+1))
      H(8)=0
      DO 395 J=1,IKK
      KS=IKK
      JJ=L14+(J-1)*4*JJJ
      JH=1
      II=0
      IXR=NOCT(L9+J)
      NIT=0
      JT=L12+5*(J-1)
      JO=L13+LZ*(J-1)
      IF(IXXX.EQ.1)GOTO 200
      NS=XXX(5+JT)
      DO 1910 IC=1, NS
 1910 IIA(100+IC)=0
      NG=NS*5
      A3=XXX(1+JT)
      B3=XXX(2+JT)
      C3=XXX(3+JT)
      D3=XXX(4+JT)
      I=0
      JD=1
      DO 121 I=1, NS
      IF(C3.EQ.0.) GOTO 121
      XE(I)=COORD(JJ+JD)
      YE(I)=COORD(JJ+JD+1)
      DI(900+I)=COORD(JJ+JD+2)
      JD=JD+4
  121 CONTINUE
C
C     THIS LOOP DETERMINES THE RELEVANT ELEMENTS AS THEY RELATE TO A
C     PARTICULAR ELEMENT.  THAT IS, EACH ELEMENT HAS ASSOCIATED WITH IT
C     THOSE OTHER ELEMENTS WHICH COULD POSSIBLY HIDE SOME PORTION
C     OF THE GIVEN ELEMENT.
C
      K=2**IGX
      K1=K
      K2=K
C
C     DO LOGARITHMIC SEARCH TO DETERMINE RELEVANT ELEMENTS.
C
      S=-1
      DO 131 I=1,IGX
      K=K+SIGN(RRX(I),S)
      IF(K.GT.IKK)K=IKK
      S=TGI(L6+J)-TGMT(L10+K)
      S1=TGI(L6+J)-TGMT(L10+K-1)
      IF(S*S1.LE.0.)GOTO 132
  131 CONTINUE
      K=IKK
  132 CONTINUE
      S=-1
      DO 133 I=1,IGX
      K1=K1+SIGN(RRX(I),S)
      IF(K1.GT.IKK)K1=IKK
      S=RV(L7+J)-YMIN(L1+K1)
      S1=RV(L7+J)-YMIN(L1+K1-1)
      IF(S*S1.LE.0.)GOTO 134
  133 CONTINUE
      K1=IKK
  134 CONTINUE
      S=-1
      DO 135 I=1,IGX
      K2=K2+SIGN(RRX(I),S)
      IF(K2.LE.1)K2=2
      IF(K2.GT.IKK)K2=IKK
      S=ZMI(L3+J)-ZMIN(L01+K2)
      S1=ZMI(L3+J)-ZMIN(L01+K2-1)
      IF(S*S1.LE.0.)GOTO 136
  135 CONTINUE
      K2=1
  136 CONTINUE
      I1=IKK-K2+1
C
C     RETRIEVE THE RELEVANT ELEMENTS DETERMINED FROM SCHEME 1.
C
      IF(XIND(L16+J).EQ.0)GOTO 1270
      IR=XIND(L16+J)
      VX=XIND(L16+J)
      T=ALOG(VX)
      IF(XIND(L16+J).LE.LC)GOTO 1800
      E=LC
      LG=XIND(L16+J)/LC
      MU=MOD(IR,LC)
      UX=LG+(MU/E)
      T=ALOG(UX)+ALOG(E)
 1800 CONTINUE
      IXT=0
      IEXP=(T/ALOG(XMD))+1
      DO 8004 L=1,IEXP
      IV=IR/(MND**(IEXP-L))
      IR=IR-IV*(MND**(IEXP-L))
      IV=IV+1
      IF(ICCT(IV-1).EQ.0)GOTO 4000
      IF(ICCT(IV-1).GT.0)GOTO 4001
      GOTO 1270
 4001 CONTINUE
      KE=ICCT(IV-1)
      IL=0
      JTT=(IV-2)*IGY+L17
      DO 4003 I=1,KE
      KV=NEH(I+JTT)
      IF(KEEP(L18+KV).EQ.J)GOTO 4003
      IL=IL+1
      NNO(L4+IXT+IL)=KV
      KEEP(L18+KV)=J
 4003 CONTINUE
      IXT=IXT+IL
 4000 CONTINUE
      IX=IBEG(IV)
      IX1=IEND(IV)
      DO 1170 I=IX,IX1
 1170 NNO(L4+IXT+I-IX+1)=IND(L15+I)
      IXT=IXT+IX1-IX+1
 8004 CONTINUE
      KS=IXT
 1270 CONTINUE
      IM=MIN0(I1,K,K1)
C
C     PICK MINIMUM COUNT FROM BOTH SCHEMES.
C
      IF(KS.LT.IM)GOTO 129
      IF(IM.EQ.I1)GOTO 1000
      IF(IM.EQ.K)GOTO 1001
      IF(IM.EQ.K1)GOTO 1002
 1000 CONTINUE
      KS=I1
      DO 1003 I=1,KS
 1003 NNO(L4+I)=IN2(L00+IKK-I+1)
      GOTO 129
 1001 CONTINUE
      KS=K
      DO 1004 I=1,KS
 1004 NNO(L4+I)=IN(L11+I)
      GOTO 129
 1002 CONTINUE
      KS=K1
      DO 1006 I=1,KS
 1006 NNO(L4+I)=IN1(L0+I)
  129 CONTINUE
      DO 170 I=1,KS
      IT=0
      MG=0
      JB=NNO(L4+I)
      IF(J.EQ.JB)GOTO 170
      JK=L13+LZ*(JB-1)
      JS=L12+5*(JB-1)
      IF((TGM(L5+J).GE.TGI(L6+JB)).OR.(TGI(L6+J).LE.TGM(L5+JB)))GOTO 170
      IF((RV(L7+J).LE.RVI(L8+JB)).OR.(RVI(L8+J).GE.RV(L7+JB)))GOTO 170
      IF(ZMI(L3+J).GE.ZM(L2+JB))GOTO 170
      NV=XXX(5+JS)
      IF(XXX(3+JS).EQ.0.) GOTO 170
      IF(XXX(3+JT).EQ.0.)GOTO 165
      NB=5*NV
C
C     TEST TO SEE IF ALL VERTICES LIE EITHER BEHIND OR IN FRONT OF
C     THE GIVEN POLYGON.
C
      JD=1
      JI=L14+(JB-1)*4*JJJ
      DO 145 M=1, NV
      A=COORD(JI+JD)
      B=COORD(JI+JD+1)
      ZI(900+M,1)=COORD(JI+JD+2)
      JD=JD+4
      XU(M)=A
      YU(M)=B
      ZS=ZI(900+M,1)
      VX=XXX(4+JT)
      VX1=XXX(2+JT)*B
      VX2=XXX(1+JT)*A
      ZS1=-(VX+VX1+VX2)/XXX(3+JT)
      IF(ABS(ZS-ZS1).LT.HXX)GOTO 145
      IT=IT+1
      ICOUNT(IT)=0
      IF(ZS.GT.ZS1)ICOUNT(IT)=1
  145 CONTINUE
      DO 1630 L=1, NS
      DO 1580 IX=1, NV
      T=ABS(YU(IX)-YE(L))
      S=ABS(XU(IX)-XE(L))
      IF((T.GT.HXX).OR.(S.GT.HXX)) GOTO 1580
      IF(ABS(ZI(900+IX,1)-DI(900+L)).GT.HXX)GOTO 1580
      IF(MG.EQ.2) GOTO 1640
      MG=MG+1
      I2(MG)=L
      I3(MG)=IX
      GOTO 1630
 1580 CONTINUE
 1630 CONTINUE
 1640 CONTINUE
      IF(MG.NE.2) GOTO 1650
      IF(J.LT.JB) GOTO 1650
      IR=IABS(I3(2)-I3(1))
      IR1=IABS(I2(2)-I2(1))
      IF((IR.NE.1).AND.(IR.NE.NV-1)) GOTO 1650
      IF((IR1.NE.1).AND.(IR1.NE.NS-1)) GOTO 1650
      IX=MAX0(I2(1),I2(2))
      IF(IABS(I2(1)-I2(2)).EQ.1)IX=IX-1
      IIA(100+IX)=1
 1650 CONTINUE
C
C     TESTS FOR SEMI-RELEVANT PLANES.  THAT IS,NEGATIVE INDEXES
C     INDICATE ELEMENT IS TO BE USED FOR VISIBILITY TEST, BUT NOT FOR
C     INTERSECTION LINE DETERMINATION.
C
      IF(IT.EQ.0)GOTO 170
      L=0
      DO 150 M=1,IT
  150 L=L+ICOUNT(M)
      IF(L.EQ.0)GOTO 170
      IF(L.EQ.IT)JB=-JB
      IF(II.NE.0)GOTO 165
C
C     INTERROGATE THE RELATIONSHIP OF THE CANDIDATE POLYGON TO THE
C     GIVEN POLYGON BY DETERMINING IF THE PROJECTION OF ONE POLYGON
C     CAN BE SEPARATED BY AN EDGE FROM THE OTHER'S PROJECTION
C
      C3=XXX(3+JT)
      C4=XXX(3+JS)
      SD=0
      I3(1)=JK
      I3(2)=JO
      I2(1)=NV*5
      I2(2)=NS*5
      DO 164 KU=1,2
      IS=I3(KU)
      IB=I2(KU)
      DO 163 L=1,IB,5
  151 CONTINUE
      IF(SD.EQ.1.)GOTO 152
      A=XXX(2+JT)*C4-XXX(2+JS)*C3
      B=XXX(1+JT)*C4-XXX(1+JS)*C3
      C=XXX(4+JT)*C4-XXX(4+JS)*C3
      GOTO 153
  152 CONTINUE
      A=CCC(L+IS)
      B=CCC(L+IS+1)
      C=CCC(L+IS+2)
  153 CONTINUE
      IF((A.EQ.0.).AND.(B.EQ.0.))GOTO 162
      IF(A.NE.0.)GOTO 154
      A=0
      C=C/B
      B=1
      GOTO 155
  154 CONTINUE
      B=B/A
      C=C/A
      A=1
  155 CONTINUE
      M=0
      R1=0
      DO 158 IX=1,NV
      M=M+1
      YG=YU(M)
      IF(A.NE.0.)GOTO 156
      DY=-C/B
      YG=XU(M)
      GOTO 157
  156 CONTINUE
      DY=-C-B*XU(M)
  157 IF(ABS(DY-YG).LT.HXX)GOTO 158
      R=YG-DY
      IF(R*R1.LT.0.)GOTO 162
      R1=R
  158 CONTINUE
      M=0
      R2=0
      DO 161 IX=1,NS
      M=M+1
      YG=YE(M)
      IF(A.NE.0.)GOTO 159
      DY=-C/B
      YG=XE(M)
      GOTO 160
  159 CONTINUE
      DY=-C-B*XE(M)
  160 IF(ABS(DY-YG).LT.HXX)GOTO 161
      R=YG-DY
      IF(R*R2.LT.0.)GOTO 162
      R2=R
  161 CONTINUE
      IF(R1*R2.LT.0.)GOTO 170
  162 CONTINUE
      IF(SD.NE.0.)GOTO 163
      SD=1
      GOTO 151
  163 CONTINUE
  164 CONTINUE
  165 CONTINUE
      II=II+1
      NNO(L4+II)=JB
  170 CONTINUE
      JS=1
      JAT=-4
      JT=L12+(J-1)*5
      NN=XXX(JT+5)
      VX=XXX(JT+4)
      VX1=XXX(2+JT)
      VX2=XXX(1+JT)
      VX3=XXX(3+JT)
      IF(IXR.LE.2)GOTO 200
      IF(II.EQ.0)GOTO 190
C
C     CALL ROUTINE WHICH SOLVES FOR THE LINES OF INTERSECTION,IF ANY,
C     OF THE JTH ELEMENT WITH OTHER ELEMENTS.
C
      CALL ISOL3D(IXR,J,XXX,CCC,II,NNO,NIT,X21,Y21,Z21,IIA,NC,ZM,ZMI,LZ)
  190 CONTINUE
  200 CONTINUE
      DO 210 JM=1,IXR
      X21(JM)=COORD(JH+JJ)
      Y21(JM)=COORD(JH+1+JJ)
      Z21(JM)=COORD(JH+2+JJ)
      IIA(JM)=COORD(JH+3+JJ)
      JH=JH+4
  210 CONTINUE
      IXR=IXR+3*NIT
      IF(II.EQ.0)GOTO 220
      IF(IXXX.NE.1)GOTO 240
  220 CONTINUE
      DO 230 JM=1,IXR
      X1(2)=X21(JM)
      Y1(2)=Y21(JM)
      IM=IIA(JM)
      CALL HPLT(X1,Y1,IJ,IM)
  230 CONTINUE
      GOTO 390
  240 CONTINUE
      JX=1
  250 CONTINUE
C
C     PLOTS IF IIA(JX+1) IS EQUAL TO 1.
C
      IF((IIA(JX).EQ.0).AND.(IIA(JX+1).EQ.0))GOTO 260
      IM=IIA(JX+1)
      X1(2)=X21(JX+1)
      Y1(2)=Y21(JX+1)
      CALL HPLT(X1,Y1,IJ,IM)
      JX=JX+2
      IF(JX.GE.IXR)GOTO 390
      GOTO 250
  260 CONTINUE
      IJ = 0
      JAT=JAT+5
      ME=0
      IF(JX.GT.NS) GOTO 255
      IF(IIA(100+JX).EQ.1) GOTO 381
  255 CONTINUE
C
C     CALL ROUTINE WHICH DETERMINES THE POINTS OF INTERSECTIONS
C     OF THE LINES OF THE JTH SET WITH THE RELEVANT LINES AND PLANES
C     OF OTHER ELEMENTS.
C
      CALL CHECK3D(XXX,CCC,NNO,J,II,NC,XI,YI,NGX,ZM,ZMI,RV,RVI,TGM,TGI,
     1	ZI,LZ)
      NG=NGX(JS)+2
      XI(1,JS)=X21(JX)
      YI(1,JS)=Y21(JX)
      ZI(1,JS)=Z21(JX)
      XI(NG,JS)=X21(JX+1)
      YI(NG,JS)=Y21(JX+1)
      ZI(NG,JS)=Z21(JX+1)
      IF(NG.LE.3)GOTO 340
C
C     THE FOLLOWING CODE SORTS THE INTERSECTION POINTS IN ASCENDING
C     ORDER OF OCCURENCE AND THEN SHRINKS THE LIST IF REDUNDANCY EXIST.
C
      NI=NG-2
      NII=NI
      DO 270 M=1,NG
      DI(M)=(XI(M,JS)-XI(1,JS))**2
      PPPPP=(YI(M,JS)-YI(1,JS))**2
      DI(M)=DI(M)+PPPPP
  270 CONTINUE
      DO 290 M=2,NI
      DO 280 MX=2,NII
      IF(DI(MX).LE.DI(MX+1))GOTO 280
      HOLD=DI(MX)
      HOLD1=XI(MX,JS)
      HOLD2=YI(MX,JS)
      HOLD3=ZI(MX,JS)
      XI(MX,JS)=XI(MX+1,JS)
      YI(MX,JS)=YI(MX+1,JS)
      ZI(MX,JS)=ZI(MX+1,JS)
      DI(MX)=DI(MX+1)
      DI(MX+1)=HOLD
      XI(MX+1,JS)=HOLD1
      YI(MX+1,JS)=HOLD2
      ZI(MX+1,JS)=HOLD3
  280 CONTINUE
      NII=NII-1
  290 CONTINUE
      LX=1
      NPX=NG
  300 NPX=NPX-1
      I=LX
      DO 320 M=I,NPX
      RX=0
      T=XI(M,JS)-XI(M+1,JS)
      T1=YI(M,JS)-YI(M+1,JS)
      T=(T**2+T1**2)**.5
      IF(T.GT.HX1)GOTO 320
      IX=M
      IX1=NPX
      DO 310 MX=IX,IX1
      XI(MX,JS)=XI(MX+1,JS)
      YI(MX,JS)=YI(MX+1,JS)
      ZI(MX,JS)=ZI(MX+1,JS)
  310 CONTINUE
      RX=1
      LX=M
      IF(LX.EQ.NPX)GOTO 330
      GOTO 300
  320 CONTINUE
  330 CONTINUE
      IF(RX.EQ.1.)NPX=NPX-1
      NG=NPX+1
  340 CONTINUE
C
C     THIS CODE DETERMINES THE STATUS(VISIBILITY) OF EVERY OTHER POINT
C     AS SUGGESTED BY THE THEOREM IN THE TECHNICAL REPORT.
C
      DO 350 L=1,NG,2
C
      OJ=XI(L,JS)
      TMJ=YI(L,JS)
      ZJ=ZI(L,JS)
      CALL STAT3D(OJ,TMJ,XXX,TGM,RV,RVI,TGI,ZM,NNO,II,H,IM,JXT,ZJ,NC,
     1	ZMI,CCC,LZ)
      DI(L)=IM
  350 CONTINUE
      DO 370 L=1,NG,2
      IF(L.EQ.NG)GOTO 370
      IF(L.EQ.NG-1)GOTO 360
      C=DI(L)+DI(L+2)
      IF(C.NE.2.)GOTO 360
      DI(L+1)=DI(L)
      GOTO 370
  360 OJ=XI(L+1,JS)
      TMJ=YI(L+1,JS)
      ZJ=ZI(L+1,JS)
      CALL STAT3D(OJ,TMJ,XXX,TGM,RV,RVI,TGI,ZM,NNO,II,H,IM,JXT,ZJ,NC,
     1	ZMI,CCC,LZ)
      DI(L+1)=IM
  370 CONTINUE
C
C     THE FOLLOWING CODE ACTUALLY PLOTS THE POINTS ON A GIVEN LINE
C     GOVERNED BY THE VALUE(IM) RETURNED BY STATUS ROUTINE.
C     1 MEANS INVISIBLE,...0 MEANS VISIBLE.
C
      DO 380 L=1,NG
      X1(2)=XI(L,JS)
      Y1(2)=YI(L,JS)
      IM=DI(L)
      CALL HPLT(X1,Y1,IJ,IM)
      IF(L.EQ.NG)GOTO 380
      C=DI(L)+DI(L+1)
      IF(C.GT.0.)GOTO 380
      OJ=(XI(L,JS)+XI(L+1,JS))/2
      TMJ=(YI(L,JS)+YI(L+1,JS))/2
      ZJ=(ZI(L,JS)+ZI(L+1,JS))/2
      CALL STAT3D(OJ,TMJ,XXX,TGM,RV,RVI,TGI,ZM,NNO,II,H,IM,JXT,ZJ,NC,
     1	ZMI,CCC,LZ)
      X1(2)=OJ
      Y1(2)=TMJ
      CALL HPLT(X1,Y1,IJ,IM)
  380 CONTINUE
  381 CONTINUE
      JX=JX+1
      GOTO 250
  390 CONTINUE
C
C     DECREMENTS THE COUNT OF THE NUMBER OF LINES IN THE JTH SET
C     SINCE THE LINES OF INTERSECTIONS WERE ADDED TO THIS ELEMENT
C     BY THE ROUTINE SOLVE.
C
      XXX(5+JT)=XXX(5+JT)-NIT
  395 CONTINUE
  400 CONTINUE
C
C     RESET VALUE FOR MAXIMUM NUMBER OF EDGES IF ARGUMENT IS
C     COMPLETED.
C
      IF (VP.GT.0.) LZ=LZ/5
      RETURN
      END
C
      SUBROUTINE STAT3D(OJ,TMJ,XXX,TGM,RV,RVI,
     1TGI,ZM,NNO,II,H,IM,JXT,ZJ,NC,ZMI,CCC,LZ)
C
C     THIS ROUTINE DETERMINES THE VISIBILITY OF AN ARBITRARY POINT
C     BY DRAWING A LINE FROM THE POINT IN QUESTION TO INFINITY AND
C     COUNTING THE NUMBER OF TIMES IT CROSSES THE BOUNDARIES OF A
C     RELEVANT ELEMENT.
C
      DIMENSION CCC(1),XXX(1)
      DIMENSION ZMI(1),TGM(1),RV(1),RVI(1),
     1	TGI(1),ZM(1),NNO(1),H(1)
C
      COMMON/GO3/L0,L1,L00,L01,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13
C
      GGK=.015
      EI=0
      IM=0
   10 CONTINUE
      IF(EI.GE.1.) GOTO 70
      EI=EI+.2
      D=EI*OJ-TMJ
      DO 60 JO=1,II
      GGK=.015
      I=0
      JG=NNO(L4+JO)
      JS=L13+(JG-1)*LZ
      JT=L12+(JG-1)*5
C
C     PRELIMINARY CHECK TO SEE IF THE POINT IS OUTSIDE THE BOUNDARY
C     BOXES IN THE X,Y,Z DIMENSIONS.
C
      IF((TMJ.GE.RV(L7+JG)).OR.(TMJ.LE.RVI(L8+JG)))GOTO 60
      IF((OJ.GE.TGI(L6+JG)).OR.(OJ.LE.TGM(L5+JG)))GOTO 60
      IF(ZJ.GE.ZM(L2+JG))GOTO 60
      VX=XXX(4+JT)
      VX1=XXX(2+JT)*TMJ
      VX2=XXX(1+JT)*OJ
      ZS=-(VX+VX1+VX2)/XXX(3+JT)
      IF(ABS(ZJ-ZS).LT.GGK)GOTO 60
      IF(ZJ.GE.ZS)GOTO 60
      NS=XXX(5+JT)
      IB=NS*5
      DO 20 J=1,IB,5
      GGK=.015
      IF(ABS(CCC(J+1+JS)).GE.100.)GGK=ALOG10(ABS(CCC(J+1+JS)))
      VE=OJ
      IF(CCC(J+JS).EQ.0.)VE=TMJ
      S=VE-CCC(J+3+JS)
      S1=VE-CCC(J+4+JS)
      YG=TMJ
      IF(CCC(J+JS).NE.0.)GOTO 15
      DY=-CCC(J+2+JS)/CCC(J+1+JS)
      YG=OJ
      GOTO 16
   15 CONTINUE
      DY=-CCC(J+2+JS)-CCC(J+1+JS)*OJ
   16 CONTINUE
      IF((ABS(YG-DY).LT.GGK).AND.(S*S1.LE.0.))GOTO 60
   20 CONTINUE
   25 CONTINUE
C
C     THE FOLLOWING CODE COUNTS THE INTERSECTIONS OF BOUNDARIES
C     OF A GIVEN ELEMENT WITH THE INFINITE LINE AND CHECKS,IF INSIDE
C     OF THE BOUNDARY, WHETHER OR NOT THE POINT IS BEHIND OR IN FRONT
C     OF THE ELEMENT.
C
      DO 40 J=1,IB,5
      T=-CCC(J+2+JS)+CCC(J+JS)*D
      R=EI*CCC(J+JS)+CCC(J+1+JS)
      IF(R.EQ.0.)GOTO 40
      T=T/R
      IF(T.LT.OJ)GOTO 40
      IF(CCC(J+JS).NE.0.)GOTO 30
      T=EI*T-D
   30 CONTINUE
      S=T-CCC(J+3+JS)
      S1=T-CCC(J+4+JS)
      IF((S.EQ.0.).OR.(S1.EQ.0.))GOTO 10
      IF(S*S1.GE.0.)GOTO 40
      I=I+1
   40 CONTINUE
      IF(MOD(I,2).EQ.0)GOTO 60
   50 CONTINUE
      IM=1
      GOTO 70
   60 CONTINUE
      IM=0
   70 CONTINUE
      RETURN
      END
C
      SUBROUTINE HPLT(X1,Y1,IJ,IM)
C
C     PLOTS POINTS GOVERNED BY THE VALUE OF IM.
C
C        NOTE THAT CALL PLOT(X,Y,2) MEANS MOVE PEN FROM THE CURRENT
C     POSITION TO THE POINT,(X,Y),WITH THE PEN DOWN.
C
C     CALL PLOT(X,Y,3) MEANS MOVE THE PEN FROM THE CURRENT POSITION
C     TO THE POINT,(X,Y), WITH THE PEN UP.
C
      DIMENSION X1(1),Y1(1)
      IF(IM.EQ.1)GOTO 20
      IF(IJ.EQ.0)GOTO 10
      CALL PLOT((X1(2)-X1(3))/X1(4),(Y1(2)-Y1(3))/Y1(4),2)
      GOTO 30
   10 CONTINUE
      CALL PLOT((X1(2)-X1(3))/X1(4),(Y1(2)-Y1(3))/Y1(4),3)
      IJ=1
      GOTO 30
   20 CONTINUE
      IJ=0
   30 CONTINUE
      RETURN
      END
C
      SUBROUTINE CHECK3D(XXX,CCC,NNO,J,II,NC,XI,YI,
     1NGX,ZM,ZMI,RV,RVI,TGM,TGI,ZI,LZ)
C
C     THIS SUBROUTINE SOLVES FOR THE POINTS OF INTERSECTION ON THE
C     LINES OF THE JTH ELEMENT WITH OTHER LINES AND PLANES(RELEVANT)
C
      DIMENSION CCC(1),XXX(1)
      DIMENSION RV(1),RVI(1),TGM(1),TGI(1),ZM(1),ZMI(1),
     1	NNO(1),NGX(15),XI(1000,1),YI(1000,1),ZI(1000,1)
C
      COMMON/DAVE/XCC(800)
      COMMON/HEDG/JS,M,JT,VX,VX1,VX2,VX3,NN
      COMMON/GO3/L0,L1,L00,L01,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13
C
      JM=1
      EEX=.015
      EXP=.005
      NGX(1)=0
      IF(II.EQ.0)GOTO 190
      IF(NN.EQ.1.) GOTO 5
      IF(VX3.NE.0.) GOTO 5
      Z=XXX(JT+2)
      B=XXX(JT+1)
      C=XXX(JT+4)
      Z1=XCC(JS)
      Z2=XCC(JS+1)
      IF(A.EQ.0.) GOTO 1
      Y1=-XCC(JS+3)*B-C
      Y2=-XCC(JS+4)*B-C
      X1=XCC(JS+3)
      X2=XCC(JS+4)
      GOTO 50
    1 CONTINUE
      Y1=XCC(JS+3)
      Y2=XCC(JS+4)
      X1=-C
      X2=X1
      GOTO 50
    5 CONTINUE
      A=XCC(JS)
      B=XCC(JS+1)
      C=XCC(JS+2)
      IF(A.EQ.0.)GOTO 20
      Y1=-XCC(JS+3)*XCC(JS+1)-XCC(JS+2)
      Y2=-XCC(JS+4)*XCC(JS+1)-XCC(JS+2)
      X1=XCC(JS+3)
      X2=XCC(JS+4)
      GOTO 30
   20 CONTINUE
      Y1=XCC(JS+3)
      Y2=XCC(JS+4)
      X1=-XCC(JS+2)
      X2=X1
   30 CONTINUE
      IF(NN.NE.1)GOTO 40
      Z1=XXX(1+JT)
      Z2=XXX(2+JT)
      GOTO 50
   40 CONTINUE
      Z1=-(VX+VX1*Y1+VX2*X1)/VX3
      Z2=-(VX+VX1*Y2+VX2*X2)/VX3
   50 CONTINUE
      AL=X2-X1
      BL=Y2-Y1
      CL=Z2-Z1
      EG=AMIN1(Z1,Z2)
      EGX=AMAX1(X1,X2)
      EGX1=AMIN1(X1,X2)
      EGY=AMAX1(Y1,Y2)
      EGY1=AMIN1(Y1,Y2)
C
C     THIS CODE DETERMINES THE POINTS OF INTERSECTIONS ON THE LINES OF
C     JTH ELEMENT RESULTING FROM THE INTERSECTION OF THE PLANES WITH
C     THESE LINES.
C
      DO 170 JR=1,II
      LG=NNO(L4+JR)
      NNO(L4+JR)=IABS(NNO(JR+L4))
      LE=NNO(L4+JR)
      JE=L13+LZ*(LE-1)
      JU=L12+5*(LE-1)
      NK=XXX(5+JU)
      JV=1
      AC=XXX(1+JU)
      BC=XXX(2+JU)
      CC=XXX(3+JU)
      D=XXX(4+JU)
      IF(EGX.LT.TGM(L5+LE))GOTO 170
      IF(EGX1.GT.TGI(L6+LE))GOTO 170
      IF(EGY.LT.RVI(L8+LE))GOTO 170
      IF(EGY1.GT.RV(L7+LE))GOTO 170
      IF(EG.GT.ZM(L2+LE))GOTO 170
      IF(LG.LT.0)GOTO 80
      IF((AL.EQ.0.).AND.(BL.EQ.0.))GOTO 80
      IF(AL.EQ.0.)GOTO 60
      XP=((BC*BL)/AL)*X1+(CC*CL/AL)*X1-D
      XP=XP-BC*Y1-CC*Z1
      VU=AC+(BC*BL/AL)+(CC*CL/AL)
      IF(VU.EQ.0.)GOTO 80
      XP=XP/VU
      T=(XP-X1)/AL
      YP=T*BL+Y1
      ZP=T*CL+Z1
      GOTO 70
   60 CONTINUE
      YP=(AC*AL/BL)*Y1+(CC*CL/BL)*Y1-D
      YP=YP-CC*Z1-AC*X1
      VU=BC+(AC*AL/BL)+(CC*CL/BL)
      IF(VU.EQ.0.)GOTO 80
      YP=YP/VU
      T=(YP-Y1)/BL
      XP=T*AL+X1
      ZP=T*CL+Z1
   70 CONTINUE
      S=ZP-ZM(L2+LE)
      S1=ZP-ZMI(L3+LE)
      IF((ABS(S).LT.EEX).OR.(ABS(S1).LT.EEX))GOTO 56
      IF(S*S1.GT.0.)GOTO 80
   56 CONTINUE
      S=XP-TGM(L5+LE)
      S1=XP-TGI(L6+LE)
      IF(S*S1.GT.0.)GOTO 80
      S=YP-RV(L7+LE)
      S1=YP-RVI(L8+LE)
      IF(S*S1.GT.0.)GOTO 80
      T=XP
      IF(A.EQ.0.)T=YP
      S=T-XCC(JS+3)
      S1=T-XCC(JS+4)
      IF(S*S1.GE.0.)GOTO 80
      M=M+1
C
C     STORES INTERSECTIONS.
C
      XI(M+1,JM)=XP
      YI(M+1,JM)=YP
      ZI(M+1,JM)=ZP
   80 CONTINUE
C
C     THIS CODE DETERMINES INTERSECTION POINTS OF LINES WITH LINES.
C
      DO 160 JC=1,NK
      B1=CCC(JV+1+JE)
      A1=CCC(JV+JE)
      C1=CCC(JV+2+JE)
      T=A1*B-B1*A
      IF(T.EQ.0.)GOTO 160
      XO=(C1*A-C*A1)/T
      IF((ABS(B).LE.50.).AND.(A.NE.0.)) GOTO 90
      YO=-C1-B1*XO
      GOTO 100
   90 CONTINUE
      YO=-C-B*XO
  100 CONTINUE
      T=XO
      IF(A.EQ.0.)T=YO
      S=T-XCC(JS+3)
      S1=T-XCC(JS+4)
      IF(S*S1.GE.0.)GOTO 160
      T=XO
      IF(A1.EQ.0.)T=YO
      S1=T-CCC(JV+4+JE)
      S=T-CCC(JV+3+JE)
      IF((ABS(S).LE.EEX).OR.(ABS(S1).LE.EEX)) GOTO 110
      IF(S*S1.GT.0.)GOTO 160
  110 CONTINUE
      IF(CC.EQ.0.)GOTO 160
      ZX=-(AC*XO+BC*YO+D)/CC
      IF((NN.NE.1.).AND.(VX3.NE.0.)) GOTO 130
      TSZ=Z2-Z1
      TSX=X2-X1
      VT=XO-X1
      IF(TSX.NE.0.)GOTO 120
      VT=YO-Y1
      TSX=Y2-Y1
  120 CONTINUE
      ZX1=(TSZ/TSX)*VT+Z1
      GOTO 140
  130 CONTINUE
      ZX1=-(VX+VX1*YO+VX2*XO)/VX3
  140 CONTINUE
      IF(ABS(ZX-ZX1).LE.EXP) GOTO 150
      IF(ZX1.GT.ZX)GOTO 160
  150 CONTINUE
      M=M+1
C
C     STORES INTERSECTIONS.
C
      XI(M+1,JM)=XO
      YI(M+1,JM)=YO
      ZI(M+1,JM)=ZX1
  160 JV=JV+5
  170 CONTINUE
      NGX(JM)=M
  190 RETURN
      END
C
      SUBROUTINE VSRTR(A,LA,IR)
C
C   IMSL ROUTINE NAME   - VSRTR
C
C-----------------------------------------------------------------------
C
C   COMPUTER            - CDC/SINGLE
C
C   LATEST REVISION     - JANUARY 1, 1978
C
C   PURPOSE             - SORTING OF ARRAYS BY ALGEBRAIC VALUE -
C                           PERMUTATIONS RETURNED
C
C   USAGE               - CALL VSRTR (A,LA,IR)
C
C   ARGUMENTS    A      - ON INPUT, A CONTAINS THE ARRAY TO BE SORTED.
C                         ON OUTPUT, A CONTAINS THE SORTED ARRAY.
C                LA     - INPUT VARIABLE CONTAINING THE NUMBER OF
C                           ELEMENTS IN THE ARRAY TO BE SORTED.
C                IR     - VECTOR OF LENGTH LA.
C                         ON INPUT, IR CONTAINS THE INTEGER VALUES
C                           1,2,...,LA. SEE REMARKS.
C                         ON OUTPUT, IR CONTAINS A RECORD OF THE
C                           PERMUTATIONS MADE ON THE VECTOR A.
C
C   PRECISION/HARDWARE  - SINGLE/ALL
C
C   REQD. IMSL ROUTINES - NONE REQUIRED
C
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP
C
C   REMARKS      THE VECTOR IR MUST BE INITIALIZED BEFORE ENTERING
C                VSRTR.  ORDINARILY, IR(1)=1, IR(2)=2, ...,
C                IR(LA)=LA.  FOR WIDER APPLICABILITY, ANY INTEGER
C                THAT IS TO BE ASSOCIATED WITH A(I) FOR I=1,2,...,LA
C                MAY BE ENTERED INTO IR(I).
C
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.
C
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN
C                           APPLIED TO THIS CODE.  NO OTHER WARRANTY,
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      DIMENSION A(1),IR(1)
C                                  SPECIFICATIONS FOR ARGUMENTS
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER            IU(21),IL(21),I,M,J,K,IJ,IT,L,ITT
      REAL               T,TT,R
C                                  FIRST EXECUTABLE STATEMENT
      IF (LA.LE.0) RETURN
      M = 1
      I = 1
      J = LA
      R = .375
    5 IF (I.EQ.J) GOTO 45
      IF (R.GT..5898437) GOTO 10
      R = R+3.90625E-2
      GOTO 15
   10 R = R-.21875
   15 K = I
C                                  SELECT A CENTRAL ELEMENT OF THE
C                                  ARRAY AND SAVE IT IN LOCATION T
      IJ = I+(J-I)*R
      T = A(IJ)
      IT = IR(IJ)
C                                  IF FIRST ELEMENT OF ARRAY IS GREATER
C                                  THAN T, INTERCHANGE WITH T
      IF (A(I).LE.T) GOTO 20
      A(IJ) = A(I)
      A(I) = T
      T = A(IJ)
      IR(IJ) = IR(I)
      IR(I) = IT
      IT = IR(IJ)
   20 L = J
C                                  IF LAST ELEMENT OF ARRAY IS LESS THAN
C                                  T, INTERCHANGE WITH T
      IF (A(J).GE.T) GOTO 30
      A(IJ) = A(J)
      A(J) = T
      T = A(IJ)
      IR(IJ) = IR(J)
      IR(J) = IT
      IT = IR(IJ)
C                                  IF FIRST ELEMENT OF ARRAY IS GREATER
C                                  THAN T, INTERCHANGE WITH T
      IF (A(I).LE.T) GOTO 30
      A(IJ) = A(I)
      A(I) = T
      T = A(IJ)
      IR(IJ) = IR(I)
      IR(I) = IT
      IT = IR(IJ)
      GOTO 30
   25 IF (A(L).EQ.A(K)) GOTO 30
      TT = A(L)
      A(L) = A(K)
      A(K) = TT
      ITT = IR(L)
      IR(L) = IR(K)
      IR(K) = ITT
C                                  FIND AN ELEMENT IN THE SECOND HALF OF
C                                  THE ARRAY WHICH IS SMALLER THAN T
   30 L = L-1
      IF (A(L).GT.T) GOTO 30
C                                  FIND AN ELEMENT IN THE FIRST HALF OF
C                                  THE ARRAY WHICH IS GREATER THAN T
   35 K = K+1
      IF (A(K).LT.T) GOTO 35
C                                  INTERCHANGE THESE ELEMENTS
      IF (K.LE.L) GOTO 25
C                                  SAVE UPPER AND LOWER SUBSCRIPTS OF
C                                  THE ARRAY YET TO BE SORTED
      IF (L-I.LE.J-K) GOTO 40
      IL(M) = I
      IU(M) = L
      I = K
      M = M+1
      GOTO 50
   40 IL(M) = K
      IU(M) = J
      J = L
      M = M+1
      GOTO 50
C                                  BEGIN AGAIN ON ANOTHER PORTION OF
C                                  THE UNSORTED ARRAY
   45 M = M-1
      IF (M.EQ.0) RETURN
      I = IL(M)
      J = IU(M)
   50 IF (J-I.GE.11) GOTO 15
      IF (I.EQ.1) GOTO 5
      I = I-1
   55 I = I+1
      IF (I.EQ.J) GOTO 45
      T = A(I+1)
      IT = IR(I+1)
      IF (A(I).LE.T) GOTO 55
      K = I
   60 A(K+1) = A(K)
      IR(K+1) = IR(K)
      K = K-1
      IF (T.LT.A(K)) GOTO 60
      A(K+1) = T
      IR(K+1) = IT
      GOTO 55
      END
C
      SUBROUTINE ISOL3D(IXR,J,XXX,CCC,II,NNO,NIT,
     1X21,Y21,Z21,IIA,NC,ZM,ZMI,LZ)
C
C     THIS ROUTINE SOLVES FOR THE LINES OF INTERSECTION RESULTING
C     FROM THE INTERSECTIONS OF THE JTH ELEMENT WITH THE OTHER
C     RELEVANT ELEMENTS.
C
      DIMENSION XXX(1),CCC(1),NNO(1),
     1ZM(1),ZMI(1),X21(1),Y21(1),Z21(1),IIA(1),IV(2)
      DIMENSION XA(50),YA(50),ZA(50)
C
      COMMON/DAVE/XCC(800)
      COMMON/GO3/L0,L1,L00,L01,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13
C
      ERS=.015
      ER=ERS
      EXX=.015
      EXP=.015
      JT=L12+(J-1)*5
      JB=L13+(J-1)*LZ
      IF(II.EQ.0)GOTO 80
      A3=XXX(1+JT)
      B3=XXX(2+JT)
      C3=XXX(3+JT)
      D3=XXX(4+JT)
      IF(XXX(3+JT).EQ.0.) GOTO 80
      DO 70 L=1,II
      K=NNO(L4+L)
C
C     CHECKS TO SEE IF THIS RELEVANT ELEMENT IS TO BE CONSIDERED FOR
C     INTERSECTION
C
      IF(K.GT.0)GOTO 9
      GOTO 70
    9 CONTINUE
      IF(K.LT.J)GOTO 70
      JX=L12+(K-1)*5
      IF(ZM(L2+J).LT.ZMI(L3+K))GOTO 70
      IF(ABS(XXX(3+JX)).LT.ERS)GOTO 70
      MT=0
      A4=XXX(1+JX)
      B4=XXX(2+JX)
      C4=XXX(3+JX)
      D4=XXX(4+JX)
C
C     DETERMINES THE EQUATION OF LINE OF INTERSECTION.
C
      B=A3*C4-A4*C3
      A=B3*C4-B4*C3
      C=D3*C4-D4*C3
      IF((A.EQ.0.).AND.(B.EQ.0.))GOTO 70
      IF(A.NE.0.)GOTO 10
      A=0
      C=C/B
      B=1
      GOTO 20
   10 CONTINUE
      B=B/A
      C=C/A
      A=1
   20 CONTINUE
      IV(1)=J
      IV(2)=K
      DO 60 M=1,2
      JV=1
      I=IV(M)
      JJ=L13+(I-1)*LZ
      IG=5+L12+(I-1)*5
      NK=XXX(IG)
      DO 50 IX=1,NK
      A1=CCC(JV+JJ)
      B1=CCC(JV+1+JJ)
      C1=CCC(JV+2+JJ)
C
C     CHECK TO BE SURE LINE OF INTERSECTION IS NOT BOUNDARY LINE
C     OF THE JTH SET.
C
      S=A1+B1+C1
      S1=A+B+C
      E=ABS(S-S1)
      S=A1*50+B1*50+C1
      S1=A*50+B*50+C
      F=ABS(S-S1)
      IF((F.LT.EXP).AND.(E.LT.EXP))GOTO 70
C
C     DETERMINES THE POINTS OF INTERSECTIONS OF THE LINE OF INTERSECTION
C     WITH OTHER LINES OF RELEVANT ELEMENTS.
C
      T=A1*B-B1*A
      IF(ABS(T).LT.ER)GOTO 50
      XO=(C1*A-C*A1)/T
      IF(A.NE.0.)GOTO 30
      YO=-C1-B1*XO
      GOTO 40
   30 CONTINUE
      YO=-C-B*XO
   40 CONTINUE
      T=XO
      IF(A1.EQ.0.)T=YO
      S=T-CCC(JV+4+JJ)
      S1=T-CCC(JV+3+JJ)
      IF(S*S1.GT.0.)GOTO 50
      MT=MT+1
C
C     STORE THE PTS OF INTERSECTIONS.
C
      XA(MT)=XO
      YA(MT)=YO
      ZA(MT)=-(D3+A3*XO+B3*YO)/C3
      ZT=-(D4+A4*XO+B4*YO)/C4
      IF(ABS(ZT-ZA(MT)).GT.EXX)GOTO 70
   50 JV=JV+5
   60 CONTINUE
      CALL STAT3D2(MT,NIT,IXR,X21,Y21,Z21,IIA,IV,A,B,
     1	C,J,XA,YA,ZA,CCC,XXX,NC,LZ)
   70 CONTINUE
   80 CONTINUE
      NR=5*XXX(5+JT)
      DO 90 IS=1,NR
      XCC(IS)=CCC(IS+JB)
   90 CONTINUE
      XXX(5+JT)=XXX(5+JT)+NIT
      RETURN
      END
C
      SUBROUTINE STAT3D2(MT,NIT,IXR,X21,Y21,Z21,
     1	IIA,IV,A,B,C,IK,XA,YA,ZA,CCC,XXX,NC,LZ)
C
C     THIS SUBROUTINE TAKES THE PTS OF INTERSECTION DETERMINED BY
C     SUBROUTINE SOLVE AND PICKS THE COORDINATES WITH THE MAX AND
C     MIN X COORDINATES PROVIDED THEY LIE ON THE INTERIOR/BOUNDARY
C     OF BOTH ELEMENTS.
C
      DIMENSION XXX(1),CCC(1),X21(1),
     1	Y21(1),Z21(1),IIA(1),IV(1),XA(1),YA(1),ZA(1)
C
      COMMON/DAVE/XCC(800)
      COMMON/GO3/L0,L1,L00,L01,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13
C
      EXX=.015
      NX=0
      IF(MT.EQ.0)GOTO 160
      DO 50 JX=1,MT
      EI=0
   10 EI=EI+.1
      IF(EI.GE.1.) GOTO 160
      D=EI*XA(JX)-YA(JX)
      DO 40 JO=1,2
      M=IV(JO)
      JC=L13+(M-1)*LZ
      JXC=L12+(M-1)*5
      NK=XXX(5+JXC)
      I=0
      IB=NK*5
C
C     DETERMINE IF THE PROJECTION OF THE POINT OF INTERSECTION
C     BELONGS TO THE INTERIOR OF BOTH PLANES.
C
      DO 30 J=1,IB,5
      EXX=.015
      IF(ABS(CCC(J+1+JC)).GE.100.)EXX=ALOG10(ABS(CCC(J+1+JC)))
      VE=XA(JX)
      IF(CCC(J+JC).EQ.0.)VE=YA(JX)
      S=VE-CCC(J+3+JC)
      S1=VE-CCC(J+4+JC)
      T=CCC(J+JC)*YA(JX)+CCC(J+1+JC)*XA(JX)+CCC(J+2+JC)
      IF((ABS(T).LT.EXX).AND.(S*S1.LE.0.))GOTO 40
      T=-CCC(J+2+JC)+CCC(J+JC)*D
      R=EI*CCC(J+JC)+CCC(J+1+JC)
      IF(R.EQ.0.)GOTO 30
      T=T/R
      IF(T.LT.XA(JX))GOTO 30
      IF(CCC(J+JC).NE.0.)GOTO 20
      T=EI*T-D
   20 CONTINUE
      IF((T.EQ.CCC(J+3+JC)).OR.(T.EQ.CCC(J+4+JC)))GOTO 10
      S=T-CCC(J+3+JC)
      S1=T-CCC(J+4+JC)
      IF(S*S1.GT.0.)GOTO 30
      I=I+1
   30 CONTINUE
      IF(MOD(I,2).EQ.0)GOTO 50
   40 CONTINUE
      NX=NX+1
      XA(NX)=XA(JX)
      YA(NX)=YA(JX)
      ZA(NX)=ZA(JX)
   50 CONTINUE
      IF(NX.EQ.0)GOTO 160
C
C     THIS CODE FINDS THE MAX/MIN X-COORDINATES(Y-COORDINATES) AND
C     STORES THEM. FUTHERMORE BOTH THE EQUATION OF LINE AND POINTS(2)
C     ARE TREATED LIKE ADDITIONAL EDGES. IN THIS WAY, THE ALGORITHM NEED
C     NOT BE DISTURBED. ESSENTIALLY,THEN,THIS TRICK IS TRANSPARENT TO
C     THE REST OF THE PROGRAM.
C
      AMAXX=-(10**6)
      AMINX=-AMAXX
      AMAXY=AMAXX
      AMINY=AMINX
      IS=5+(IK-1)*5+L12
      IS=XXX(IS)
      DO 110 JI=1,NX
      IF(A.EQ.0.)GOTO 80
      IF(XA(JI).GE.AMINX)GOTO 60
      AMINX=XA(JI)
      YI=YA(JI)
      ZI=ZA(JI)
   60 IF(XA(JI).LE.AMAXX)GOTO 70
      AMAXX=XA(JI)
      YII=YA(JI)
      ZII=ZA(JI)
   70 CONTINUE
      GOTO 110
   80 CONTINUE
      IF(YA(JI).GE.AMINY)GOTO 90
      AMINY=YA(JI)
      XI=XA(JI)
      ZI=ZA(JI)
   90 CONTINUE
      IF(YA(JI).LE.AMAXY)GOTO 100
      XII=XA(JI)
      AMAXY=YA(JI)
      ZII=ZA(JI)
  100 CONTINUE
  110 CONTINUE
      NIT=NIT+1
      K=5*(NIT-1+IS)+1
      XCC(K)=A
      XCC(K+1)=B
      XCC(K+2)=C
      IF(A.EQ.0.)GOTO 120
      XCC(K+3)=AMINX
      XCC(K+4)=AMAXX
      AMIN=AMINX
      AMAX=AMAXX
      YE=YII
      ZE=ZII
      GOTO 130
  120 CONTINUE
      XCC(K+3)=AMINY
      XCC(K+4)=AMAXY
      AMIN=XI
      AMAX=XII
      YI=AMINY
      YE=AMAXY
      ZE=ZII
  130 CONTINUE
      IG=IXR+NIT*3
      X21(IG-2)=AMIN
      Y21(IG-2)=YI
      Z21(IG-2)=ZI
      DO 140 JK=1,2
      IE=IG-JK+1
      X21(IE)=AMAX
      Y21(IE)=YE
      Z21(IE)=ZE
  140 CONTINUE
      DO 150 JK=1,2
      IIA(IG-JK)=0
  150 CONTINUE
      IIA(IG)=1
      TX=(AMAX-AMIN)**2
      TY=(YE-YI)**2
      DX=(TX+TY)**.5
      IF(DX.LT..001)NIT=NIT-1
  160 RETURN
      END
C
      SUBROUTINE VSRT1(A,LA,IR)
C
      INTEGER IU(21),IL(21),I,M,J,K,IJ,IT,L,ITT
      INTEGER A(1),IR(1),T,TT
C                                  FIRST EXECUTABLE STATEMENT
      IF (LA.LE.0) RETURN
      M = 1
      I = 1
      J = LA
      R = .375
    5 IF (I.EQ.J) GOTO 45
      IF (R.GT..5898437) GOTO 10
      R = R+3.90625E-2
      GOTO 15
   10 R = R-.21875
   15 K = I
C                                  SELECT A CENTRAL ELEMENT OF THE
C                                  ARRAY AND SAVE IT IN LOCATION T
      IJ = I+(J-I)*R
      T = A(IJ)
      IT = IR(IJ)
C                                  IF FIRST ELEMENT OF ARRAY IS GREATER
C                                  THAN T, INTERCHANGE WITH T
      IF (A(I).LE.T) GOTO 20
      A(IJ) = A(I)
      A(I) = T
      T = A(IJ)
      IR(IJ) = IR(I)
      IR(I) = IT
      IT = IR(IJ)
   20 L = J
C                                  IF LAST ELEMENT OF ARRAY IS LESS THAN
C                                  T, INTERCHANGE WITH T
      IF (A(J).GE.T) GOTO 30
      A(IJ) = A(J)
      A(J) = T
      T = A(IJ)
      IR(IJ) = IR(J)
      IR(J) = IT
      IT = IR(IJ)
C                                  IF FIRST ELEMENT OF ARRAY IS GREATER
C                                  THAN T, INTERCHANGE WITH T
      IF (A(I).LE.T) GOTO 30
      A(IJ) = A(I)
      A(I) = T
      T = A(IJ)
      IR(IJ) = IR(I)
      IR(I) = IT
      IT = IR(IJ)
      GOTO 30
   25 IF (A(L).EQ.A(K)) GOTO 30
      TT = A(L)
      A(L) = A(K)
      A(K) = TT
      ITT = IR(L)
      IR(L) = IR(K)
      IR(K) = ITT
C                                  FIND AN ELEMENT IN THE SECOND HALF OF
C                                  THE ARRAY WHICH IS SMALLER THAN T
   30 L = L-1
      IF (A(L).GT.T) GOTO 30
C                                  FIND AN ELEMENT IN THE FIRST HALF OF
C                                  THE ARRAY WHICH IS GREATER THAN T
   35 K = K+1
      IF (A(K).LT.T) GOTO 35
C                                  INTERCHANGE THESE ELEMENTS
      IF (K.LE.L) GOTO 25
C                                  SAVE UPPER AND LOWER SUBSCRIPTS OF
C                                  THE ARRAY YET TO BE SORTED
      IF (L-I.LE.J-K) GOTO 40
      IL(M) = I
      IU(M) = L
      I = K
      M = M+1
      GOTO 50
   40 IL(M) = K
      IU(M) = J
      J = L
      M = M+1
      GOTO 50
C                                  BEGIN AGAIN ON ANOTHER PORTION OF
C                                  THE UNSORTED ARRAY
   45 M = M-1
      IF (M.EQ.0) RETURN
      I = IL(M)
      J = IU(M)
   50 IF (J-I.GE.11) GOTO 15
      IF (I.EQ.1) GOTO 5
      I = I-1
   55 I = I+1
      IF (I.EQ.J) GOTO 45
      T = A(I+1)
      IT = IR(I+1)
      IF (A(I).LE.T) GOTO 55
      K = I
   60 A(K+1) = A(K)
      IR(K+1) = IR(K)
      K = K-1
      IF (T.LT.A(K)) GOTO 60
      A(K+1) = T
      IR(K+1) = IT
      GOTO 55
      END
